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are  taken  to  be  ai  =  (a/2)  (ev  +  ez),  a2  =  (a/2)  (ex+  ez),  and  a3  =  (a/2)  (ex +ey),  where  ex, 
ey,  and  ez  are  unit  vectors  pointing  along  the  x-,  y-,  and  z-directions  shown  in  subfigures 


(a)  and  (b)  (12).  Special  high-symmetry  points  in  this  Brillouin  zone  are  denoted  by  F, 

K,  W ,  X,  U,  and  L.  The  line  segments  connecting  points  F  and  K,  F  and  X,  and  T  and  L 
are  denoted  as  £,  A,  and  A,  respectively  (17) . 7 
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1.  Introduction 


Since  the  1950s,  strain  has  been  known  to  control  the  electronic  band  structure  of 
semiconductors  (1 ).  This  mechanism  was  used  as  a  parameter  to  design  novel  semiconductor 
material  technologies  recently  in  complementary  metal-oxide-semiconductor  (CMOS) 
technologies  to  control  the  mobility  of  electrons  in  short  field-effect  transistor  channels.  In 
doing  so,  it  enabled  Intel  Corporation  to  move  from  130-  to  90-nm  feature  sizes  in  2003  (2), 
which  effectively  enabled  the  continued  realization  of  Moore’s  Law  and  still  does  so  today. 

Modeling,  particularly  with  regard  to  the  advent  of  strained  silicon,  played  a  significant  role  in 
the  development  of  early  science  and  breakthrough  technologies.  Indeed,  the  basic  mechanisms 
of  strained  silicon,  and  their  theoretical  underpinnings,  were  known  as  early  as  the  1950s  in 
many  of  the  historic  papers  that  led  the  way  for  CMOS  and  p-type  metal-oxide-semiconductor 
(PMOS)  transistors  (/).  Computations  thus  had  a  theoretical  and  mathematical  framework  from 
which  to  draw  models  and  computable  concepts. 

Today,  the  majority  of  relevant  scientific  papers  fall  into  one  of  two  basic  problem  groups. 
Roughly  speaking,  the  first  focuses  on  the  determination  of  the  electronic  band  structure  and  the 
second  focuses  on  the  time-varying  nature  of  electron  transport.  Computational  research  can 
likewise  fit  within  this  loose  classification.  On  their  own,  each  represents  large  but  separate 
areas  of  scientific  endeavor.  But  together,  they  provide  a  basis  for  robust  design  of 
semiconductor  device  materials. 

With  greater  access  now  sought  for  more  chemical  species  in  the  periodic  table  of  elements  to  be 
considered  for  doping,  substrate  design,  or  even  the  main  semiconductor  material,  the  many 
computational  approaches  that  provide  an  early  testing  platform  for  materials  designs  demand  a 
reexamination  for  their  transferability.  Of  particular  connection  to  the  Army,  strain  is  a  potential 
issue  in  many  optoelectronic  devices  that  are  composed  of  multiple  materials,  such  as 
semiconductor  lasers  or  light-emitting  diodes  (LEDs).  Such  materials  and  devices  are  not 
uncommon  in  infrared,  ultraviolet,  and  microwave  detectors  and  sources  either  sought  by  or 
under  consideration  as  new  concepts  by  the  U.S.  Army. 

As  a  specific  motivation,  consider  the  optical  gain  medium  in  quantum  dot  lasers  composed  of  an 
array  of  indium  arsenide  (InAs)  inclusions  embedded  in  a  gallium  arsenide  (GaAs)  matrix  that 
have  relatively  recently  been  commercialized  (3-5).  The  mismatch  between  the  lattice  constants 
of  these  two  materials  is  substantial:  7%  ( 6 ).  Strain  may  affect  the  states  of  optically  active 
electrons  through  piezoelectric  effects,  especially  for  nitride  semiconductor  optoelectronics  (7, 

8).  However,  even  where  such  effects  are  absent  or  negligible,  strain  still  has  important  effects 
on  the  electronic  properties  of  a  device,  such  as  the  minimum  energy  needed  to  promote  an 
electron  from  the  ground  state  to  an  excited  state  (9, 10). 


1 


In  this  report,  we  perfonn  a  brief  survey  of  various  existing  methods  used  to  estimate  the 
energies  and  quantum-mechanical  wave  functions  of  systems  of  electrons,  how  these  methods 
take  strain  into  account,  and  where  they  intersect  with  concepts  from  solid  mechanics.  As  such, 
we  consider  only  the  first  part  of  the  two-part  problem — the  accurate  computational 
determination  of  electronic  band  structure  in  strained  materials  and  structures.  The  survey  will 
predominantly  cover  approximate  empirical  approaches.  However,  in  order  to  make  clear  what 
impact  the  approximations  in  these  approaches  have  had,  the  survey  will  include  a  brief 
discussion  of  the  computationally  expensive  quantum-mechanical  approaches  known  as  ab  initio 
methods,  which  will  be  contrasted  with  the  empirical  electronic  structure  methods.  The 
multiscale  aspects  of  these  empirical  methods  will  be  highlighted. 

The  intended  audience  is  the  community  of  researchers  with  interest  in  accurately  computing  the 
electronic  structure  properties  and,  in  particular,  the  excited  state  properties  of  materials. 
However,  considering  the  intersection  with  concepts  of  strain,  the  intended  community  would  be 
more  familiar  with  concepts  of  elasticity  and  deformability  of  solids.  Thus,  section  2  starts  with 
a  brief  review  of  electronic  structure  theory,  providing  the  mathematical  starting  point  to  which 
later  assumptions  are  ascribed  for  the  sake  of  computational  feasibility.  It  then  briefly  describes 
ab  initio  methods  and  tersely  lays  the  foundations  for  first-principles  approaches,  including 
discussion  of  spin-orbit  coupling.  Section  3  presents  an  overview  of  empirical  atomistic 
approaches  that  covers  the  empirical  pseudopotential  and  the  Slater-Koster  tight-binding 
methods.  Section  4  describes  the  so-called  continuum  methods,  specifically  the  interrelated 
k  ■  p  envelope  function,  and  effective  mass  methods.  A  potential  issue  for  computing  electronic 
structures  of  heterogeneous  domains  comprised  of  multiple  materials — the  valence  band  offset — 
is  described  in  section  5  within  the  context  of  the  empirical  approaches  covered  in  this  report.  In 
section  6,  we  show  computed  comparisons  of  the  different  methods  using  models  of  GaAs,  InAs, 
and  aluminum  arsenide  (AlAs).  Section  7  completes  the  report  with  conclusions. 


2.  Background:  Electronic  Structure  Theory 


In  principle,  the  electronic  structure  of  a  system  of  n  electrons  and  M  atomic  nuclei,  whether  it 
be  a  bulk  crystal,  molecule,  nanostructure,  etc.,  is  the  solution  of  a  many-body  quantum 
mechanical  differential  equation.  For  a  nonrelativistic  system,  this  equation  is  the  (time- 
independent)  many-electron  Schrodinger  equation  (11): 

Wmanye-^CWn)  =  ^Tot^CWn)  ,  (1) 

where  (r}n  is  shorthand  for  all  of  the  electron  coordinates  r1;  r2, ... ,  rn  and  their  spins 
s1;  s2, ... ,  sn.  ^({r}^)  is  the  kth  many-electron  wavefunction,  and  for  this  wavefunction,  E^otk 
is  the  total  energy  of  the  n  electrons  in  the  system.  Hmany  e-  is  a  differential  operator  called  the 
many-electron  Hamiltonian.  Here,  the  notation  is  such  that  Ejot  0  <  Ejot  l  <  Ejot2,  etc.,  where 
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E-rot, o  is  the  lowest  possible  energy  of  the  electrons  in  the  system,  that  is,  the  ground-state 
energy,  and  T0({r)n)  denotes  a  ground-state  many-electron  wavefunction.  (If,  for 
example,  Ejot0  —  Ejot  l,  then  ((r}n)  is  also  a  ground-state  wavefunction.)  Wavefunctions 
with  energies  higher  than  Ejot  0  correspond  to  excited  states  of  the  system.  The  physical 
meaning  of  is  such  that 


P  =  X 


Sl,S2,. 


snL  L  ■" /o„|h//c({r}n)l2d3r1d3r2  ...d3r„ 


(2) 


is  the  probability  that  electron  1  is  in  the  region  fi1?  electron  2  is  in  the  region  fl2,  etc.,  provided 
that  H'j  is  normalized  so  that  P  =  1  if  fl2, ...  ,fln->oo.  The  summation  in  equation  2  is  over 
all  possible  spin  values  of  the  electrons.  The  many-electron  wavefunction  must  satisfy  the  Pauli 
exclusion  principle,  which  states  that  two  electrons  cannot  share  the  same  state.  The  principle  is 
satisfied  if,  and  only  if,  this  wavefunction  changes  sign  when  two  electronic  coordinates  are 
exchanged,  that  is, 

vFi(r1>  Sp,  r2,  s2, ... ,  ri7  S;,  ri+1,  si+1 ... ,  rn,  sn) 

=  r2, ... ,  ri+1,  si+1,  ri;  S;, ... ,  rn,  sn)  .  (3) 

If  two  electrons  were  to  share  the  same  state,  i.e.,  if  fa,  sL)  =  (ri+1, si+1),  then  Tj  would  have  to 
equal  zero  for  equation  3  to  hold. 

If  the  nuclei  are  treated  as  classical  point  particles,  then  Hmanye-  takes  the  following  form  (11): 


H, 


many  e 


yn 

1*1=1 


rplpl 
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+  Vf 


ext  (r,,{R}M)]  +  Pinter({r}) 


(4) 


Here,  pl  =  —  iftVj  is  the  momentum  operator,  i  =  V— T,  h  —  h/2n  is  the  reduced  Planck 
constant,  V j  is  the  del  operator  defined  with  respect  to  electron  coordinate  me  is  the  mass  of 
an  electron,  and 


Fext(r>  (R)m) 


1  y  M  Zjq2 
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where  q  is  the  magnitude  of  the  charge  of  an  electron,  e0  is  the  permittivity  of  free  space,  and 
(R}m  is  shorthand  for  the  coordinates  of  all  M  nuclei  of  the  system  R1;  R2, ... ,  RM  and  their 
respective  atomic  numbers  Zv  Z2, ... ,  ZM.  pl  •  p l/2 me  is  the  kinetic  energy  operator  for  electron 
i.  The  sum  of  Yi=i  kCxt(r(-  (R)m)  and  Pinter ((r)n)  is  the  total  electronic  potential  energy. 

Vext(r,  (R}m)  is  called  an  external  potential  and  accounts  for  the  attractive  electrostatic 
interactions  between  each  electron  and  the  nuclei,  while  Pinter  (Mn)  is  an  interaction  potential 
accounting  for  the  repulsions  between  the  electrons. 
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The  presence  of  Einter({r}n)  prevents  the  n-electron  Schrodinger  equation  from  being  separated 
into  n  individual  equations,  which  makes  an  exact  solution  of  the  many-electron  Schrodinger 
equation  intractable  for  n  >  1.  Instead,  the  so-called  independent  electron  approximation  is 
adopted,  where  each  electron  is  nominally  treated  as  if  it  were  independent,  but  the  actual 
interaction  of  each  electron  with  all  other  electrons  in  the  system  is  taken  into  account  through  an 
effective  external  potential  (12).  This  leads  to  a  one-electron  Schrodinger  equation,  which  is 


Hie-4>i(r,  s ) 


P  P 


L2m, 


T  ^ext,eff(T’  S>  {R}m) 


xpi(r,s)  =  Etipifas) 


for  an  isolated  system,  or 


(7) 


#le-l/hkO,s) 
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Eud/hk(r,s) 


(8) 


for  a  system  that  is  periodic  along  one  or  more  dimensions.  Here,  Hle-  is  the  one-electron 
Hamiltonian,  r  and  s  are  the  position  and  spin  of  the  electron,  p  =  —  iftV  is  the  momentum 
operator,  p  •  p/2 me  is  the  kinetic  energy  operator,  Vext,eff(r»  s,  {R})  is  an  effective  external 
potential  operator,  and  ipi(r,  s),  EL  [or  ipik(r,  s),  and  Eik]  are  the  one-electron  wavefunction  and 
corresponding  eigenenergy  of  electron  i  (or  ik).  The  quantity  k  is  called  a  Bloch  wavevector, 
which  will  be  defined  and  discussed  later.  The  effective  external  potential  may  be  decomposed 
as 


^ext,eff(r, s'  (r)m)  =  ^extCr,  Wm)  +  t4_e(r,  s)  .  (9) 

Here,  l4-e(r< s)  is  an  operator  that  takes  into  account  the  interaction  between  an  electron  at  r  and 
all  the  other  electrons  in  the  system.  The  relationship  between  the  one-electron  wavefunction 
ipi(r,s),  or  xpik(r,s),  and  the  full  many-body  wavefunction  ^((r^)  depends  on  l4_e(r,.s),  so 
approximations  to  14-e(r< s)  lead  to  certain  approximations  to  this  relationship,  as  seen  in  section 
2. 1  on  ab  initio  methods.  The  eigenenergies  Et,  or  Eik,  can  be  viewed  as  a  set  of  energy  levels 
(11).  Examples  of  this  are  pictured  in  figure  1,  which  shows  a  schematic  of  these  levels  for  an 
isolated  atom  with  a  nuclear  charge  of  n  in  four  different  configurations:  (1)  ground  state  of  the 
atom  when  it  is  neutral,  i.e.,  the  number  of  electrons  equals  the  nuclear  charge  n;  (2)  ground  state 
of  the  atom  with  an  extra  negative  charge;  (3)  ground  state  of  the  atom  with  an  extra  positive 
charge;  and  (4)  a  neutral  excited  state  of  the  atom.  The  total  electronic  energies  of  each  state  are 
shown.  Because  of  the  Pauli  exclusion  principle,  only  two  electrons  (one  spin-up  and  one  spin- 
down)  can  occupy  each  energy  level.  If  n  is  even,  all  the  levels  up  to  EceiKn/2)  are  occupied  in 
the  ground  state  of  a  neutral  atom.  Otherwise,  one  of  the  energy  levels  is  partially  occupied, 
even  in  the  ground  state.  [The  ceil(  ) function  returns  the  smallest  integer  greater  than  or  equal 
to  its  argument.]  Here,  energy  levels  above  the  highest  ground-state  energy  level,  Ecei^n/ 2),  are 
energies  required  to  add  an  electron  to  the  system  at  a  certain  level,  so  Ece ii(n/2)+i  — 

~^Tot,0’  ECcii(n/2)+2  —  Etou  —  E-rot, 0’  etc-  Energies  equal  to  or  below  the  highest  ground-state 
energy  level  are  energies  required  to  remove  an  electron  from  the  system,  so  Eceil(n/2) 

=  E-rot  0  —  E-po^o ,  Eceil(n/2)_i  =  E-rot  0  —  Ej^,  etc.  (13).  To  a  first  approximation,  the  available 
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Figure  1.  Energy  levels  of  an  isolated  atom  with  a  nuclear  charge  of  n  in  four 
different  configurations:  (1)  ground  state  of  the  atom  when  it  is 
neutral,  i.e.,  the  number  of  electrons  equals  the  nuclear  charge  n, 

(2)  ground  state  of  the  atom  with  an  extra  negative  charge,  (3)  ground 
state  of  the  atom  with  an  extra  positive  charge,  and  (4)  a  neutral 
excited  state  of  the  atom,  where  an  electron  has  been  excited  from 
level  ceil(«/2)  to  ceil(«/2)+l.  The  total  electronic  energies  of  each 
state  are  shown.  The  ceil()  function,  which  returns  the  smallest 
integer  greater  than  or  equal  to  its  argument,  accounts  for  the  case 
where  n  is  odd.  In  the  diagram  shown,  however,  n  is  even. 

neutral  excited  states  may  be  estimated  from  these  energy  levels  that — strictly  speaking — pertain 
to  addition  and  removal  of  electrons.  For  example,  a  neutral  excited  state  may  be  treated  as  if  it 
were  due  to  an  addition  of  an  electron  at  a  higher  energy  level  followed  by  the  removal  of  an 
electron  at  a  lower  level.  Such  a  state  is  shown  in  figure  1,  where  an  electron  is  excited  to  the 
next  higher  energy  level  and  increased  in  energy  by  FCeii(n/2)+i  —  Ecei](jl/2y  An  electron  in  such 
an  excited  state  will  emit  a  photon  with  that  amount  of  energy  when  it  returns  to  the  ground  state, 
a  process  that  is  important  in  devices  such  as  lasers  and  LEDs  (14).  Naturally,  however,  this 
approximation  is  best  suited  for  the  addition  or  removal  of  electrons  in  the  system,  as  opposed  to 
elementary  excitations  such  as  bound  states  or  electron-hole  pairs  (13, 15, 16).  Additional 
energy  levels  may  be  introduced  by  the  attraction  between  an  excited  electron  and  the  vacancy  or 
“hole”  it  leaves  behind  in  its  former  energy  level.  These  new  levels  cannot  be  predicted  from  the 
addition  and  removal  energies  of  the  independent-electron  approximation.  For  simplicity,  the 
energy  levels  in  the  figure  have  been  presented  as  if  they  were  distinct,  but  this  is  often  not  the 
case.  Rather,  different  one-electron  wavefunctions,  e.g.,  xpt  and  xpi+1,  can  share  the  same  energy 
eigenvalue;  that  is,  Et  —  Et+1.  Such  energy  levels  are  called  degenerate. 

For  periodic  systems,  the  one-electron  eigenenergy  Eik  depends  on  the  quantity  k  in  equation  8, 
which  comes  from  the  Bloch  theorem.  This  theorem  states  that  if  the  potential  Fext,eff  is 
periodic,  the  one-electron  wavefunction  has  the  form  (12) 
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*/>ik(r,s)  =  eikritik(r,  s)  ,  (10) 

where  k  is  called  a  Bloch  wavevector,  and  uik(r,  s)  has  the  same  periodicity  as  Fext,eff-  The  one- 
electron  wavefunction,  then,  has  the  form  of  a  traveling  wave  propagating  in  the  direction  of  k. 
The  periodicity  of  the  potential  can  be  expressed  in  terms  of  a  lattice  vector  R  =  nqaj, 
where  m £  and  Np  are  integers  and  a,  is  a  primitive  lattice  vector.  Accordingly, 

Text,eff(T>  Text,eff(r  +  R,  S,  {R}m)  , 

and  (11) 

uik (r,s)  =  uik(r  +  R,s) . 

If  the  system  is  periodic  in  all  three  directions,  then  Np  —  3  and  the  three  lattice  vectors  form  a 
unit  cell,  as  shown  in  figure  2.  The  Bloch  wave  vector  can  then  point  along  any  direction  in 
space.  If  the  periodicity  is  confined  to  a  plane,  then  Np  =  2  and  the  wave  vector  k  is  also 
confined  within  the  plane  defined  by  the  two  lattice  vectors  a1  and  a2.  If  there  is  only 
periodicity  along  one  dimension,  then  Np  —  1  and  k  is  parallel  or  antiparallel  to  a1  (10).  For  any 
value  of  Alp,  a  reciprocal  space  can  be  defined  such  that  its  primitive  vectors  satisfy  the 
relationship  bj  •  a;  =  2n8Lj,  where  is  the  Kronecker  delta.  A  general  vector  in  reciprocal 
space  is  q  =  q  bj,  where  q  may  be  any  real  number;  a  general  reciprocal  lattice  vector, 
then,  is  G  =  Yn=ini  bj,  where  nj  is  an  integer.  In  general,  elRG  =  1.  The  Bloch  wave  vector  k 
exists  within  the  subset  of  reciprocal  space  called  either  the  first  Brillouin  zone  or  simply  the 
Brillouin  zone,  which  consists  of  the  points  in  reciprocal  space  closer  to  G  =  0  than  to  any  other 
reciprocal  lattice  point  (12). 


Figure  2.  The  three  primitive  lattice  vectors  of  3-D 

crystalline  unit  cell,  which  is  shown  in  dashed 
lines  as  a  parallelepiped.  The  origin  of  the 
primitive  lattice  vectors  is  shown  as  located  at 
a  corner  of  a  unit  cell  of  a  crystal,  but  the 
origin  is  arbitrary  and  can  be  taken  to  be,  for 
example,  the  center  of  a  unit  cell. 
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Figure  3d  shows  the  Brillouin  zone  of  a  bulk  crystal  with  a  zincblende  or  diamond-type  structure. 
These  two  structures  are  common  in  semiconductors  such  as  GaAs  and  silicon  (Si).  The 
conventional  cubic  unit  cells  of  such  crystal  structures  are  shown  in  figures  3a  and  b,  while  the 
primitive  unit  cells  (i.e.,  the  smallest  possible  unit  cells  needed  to  specify  the  crystal  structure) 
are  shown  in  figure  3c.  Certain  high-symmetry  points  in  the  Brillouin  zone  are  given  special 
labels.  For  example,  the  center  of  the  zone,  where  k  =  0,  is  denoted  a  .  The  labels  for  other 
high-symmetry  points  are  shown  in  figure  3d.  Typically,  Elk  is  plotted  for  values  of  k  that  trace 
a  path  connecting  several  of  these  high-symmetry  points  in  the  Brillouin  zone  ( 1 7),  and  such  a 
plot  is  shown  in  figure  4,  where  the  path  traced  is  from  L  to  o  to  X  to  K  and  back  to  o  .  These 
plots  are  diagrams  of  the  band  structure  of  the  crystal. 


(d) 


Figure  3.  (a)  Conventional  unit  cell  of  a  material  with  a  zincblende  crystal  structure  and  lattice  constant  a. 

(b)  Conventional  unit  cell  of  a  material  with  a  diamond  crystal  structure  and  a  lattice  constant  a.  The 
only  difference  between  the  zincblende  and  diamond  structures  is  that  the  atoms  in  the  latter  are  all  of 
the  same  type,  (c)  The  primitive  unit  cells  for  zincblende-  and  diamond-type  crystals,  (d)  Brillouin 
zone  corresponding  to  the  primitive  cells  shown  in  subfigure  (c).  The  primitive  lattice  vectors  are  taken 
to  be  ai  =  (a/2)  (ev  +  ez),  a2  =  (a/2)  (ez+  ez),  and  a3  =  (a/2)  (e,.+ey),  where  e  ,  ey,  and  ezare  unit  vectors 
pointing  along  the  x-,  y-,  and  z-directions  shown  in  subfigures  (a)  and  (b)  (12).  Special  high-symmetry 
points  in  this  Brillouin  zone  are  denoted  by  a  ,  K,  W.  X,  U,  and  L.  The  line  segments  connecting  points 
a  and  K,  a  and  X,  and  a  and  L  are  denoted  as  G  ,  r  ,  and  v  ,  respectively  (77). 


7 


(a)  (b) 


Figure  4.  Band  structures  of  (a)  GaAs  and  (b)  Si,  as  estimated  by  the  method  of  Vogl  et  al.  ( 18).  In 

each  diagram,  the  locations  of  the  valence  band  maximum  and  conduction  band  minimum  are 
shown.  The  zero-energy  datum  is  taken  to  be  the  valence  band  maximum. 

As  mentioned  before,  the  energy  levels  of  an  isolated  n-electron  atom  in  its  ground  state  are 
occupied  only  up  to  a  certain  maximum  level  Ecei^n/2).  Similarly,  in  a  crystal  in  its  ground  state, 
only  bands  with  energies  below  a  certain  maximum  level  are  occupied.  For  semiconductors  and 
insulators,  this  maximum  level  is  the  valence  band  maximum.  Figure  4  shows  the  valence  band 
maxima  for  Si  and  GaAs.  For  an  electron  in  these  materials  to  be  excited  to  the  next  highest 
band,  called  the  conduction  band,  its  energy  must  become  at  least  equal  to  the  conduction  band 
minimum,  which  is  also  shown  in  figure  4.  The  difference  between  the  valence  band  maximum 
and  conduction  band  minimum  is  called  the  band  gap  energy  and  is  denoted  here  as  Eg;  it 
governs  the  frequency  of  the  photon  emitted  by  an  excited  electron  as  it  returns  to  the  ground 
state.  In  the  band  structure  for  GaAs  shown  in  figure  4a,  the  conduction  band  minimum  and 
valence  band  maximum  occur  for  the  same  Bloch  wave  vector  k  =  0,  so  GaAs  is  called  a  direct 
band  gap  material.  All  that  is  needed,  then,  for  an  electron  with  such  a  wave  vector  to  be 
promoted  from  the  valence  to  the  conduction  band  is  to  absorb  an  photon  with  energy  Eg.  The 
band  structure  for  Si  shown  in  figure  4b  shows  the  conduction  band  minimum  and  valence  band 
maximum  occurring  at  different  wave  vector  values.  In  such  a  case,  to  promote  an  electron  from 
the  valence  band  to  the  conduction  band  it  is  not  enough  for  an  electron  to  simply  absorb  a 
photon  with  an  energy  equal  to  Eg .  The  electron  must  also  have  momentum  imparted  to  it  by 
thennal  vibrations  of  the  crystal  nuclei  in  order  to  change  its  wave  vector  to  that  of  the 
conduction  band  minimum.  For  this  reason,  Si  is  called  an  indirect  band  gap  material  (12). 
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Metals,  unlike  semiconductors  and  insulators,  have  no  band  gap  at  all,  as  seen  in  the  band 
structure  of  silver  (Ag)  in  figure  5,  which  is  why  metals  readily  conduct.  This  illustrates  the 
importance  of  band  structure  in  determining  electronic  properties. 


k 


Figure  5.  Band  structure  of  Ag,  as  determined  via  the  Naval  Research  Laboratory  (NRL) 
tight-binding  code  (19,  20).  Since  the  lattice  vectors  of  the  primitive  cells  of  the 
diamond  and  zincblende  structures  are  the  same  as  that  of  the  primitive  cell  of 
the  crystal  structure  of  Ag  (face-centered  cubic,  or  fee),  the  special  points  along 
the  horizontal  axis  of  this  diagram  are  the  same  as  those  in  figure  3d.  The 
primitive  cell  of  an  fee  lattice  has  only  one  unit  cell. 

To  summarize,  while  the  electronic  structure  of  a  system  is  in  principle  the  solution  of  an 
intractable  many-body  equation,  in  practice  it  can  be  characterized  by  a  picture  where  the 
electrons  are  nominally  treated  as  independent  and  each  have  their  own  eigenenergies  and  one- 
electron  wavefunctions.  For  nonperiodic  systems,  these  eigenenergies  are  discrete  values  Et,  but 
periodicity  in  one  or  more  directions  introduces  a  wavevector  k  on  which  the  eigenenergies,  now 
denoted  as  Eik,  also  depend.  Examination  of  these  eigenenergies  can  shed  light  on  various 
physical  properties  of  the  system,  as  shown  for  relatively  simple  examples  such  as  isolated  atoms 
and  bulk  crystals.  How  one  may  obtain  an  effective  one-electron  equation  from  the  many- 
electron  equation,  and  the  relationship  of  the  one-electron  wavefunctions  to  the  many-body 
wavefunction,  will  be  discussed  in  the  following  section  on  ab  initio  methods. 

2,1  Ab  Initio  Methods 

Ab  initio  methods  amount  to  ways  of  either  approximately  solving  the  many-electron 
Schrodinger  equation  or  detennining  quantities  of  interest,  such  as  £'Tot,o<  without  directly 
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solving  the  many-electron  equation  itself.  Often,  these  methods  involve  solving  an  effective 
one-electron  Schrodinger  equation.  Some  methods  are  characterized  by  the  approximations 
involved  in  formulating  the  operator  Vc_c  that  take  into  account  the  interaction  between 
electrons.  Other  formulations  may  be  exact  in  principle  but  still  require  approximations  to  obtain 
solutions  in  practice.  For  the  sake  of  simplicity,  dependencies  on  the  wave  vector  k  or  spin  are 
suppressed  in  the  following  briefly  outlined  methods. 


One  of  the  earliest  ab  initio  approaches,  the  Hartree  approximation,  makes  calculations  tractable 
through  the  assumption  that  the  electrons  of  the  system  can  be  treated  as  if  they  were  smeared 
out  into  a  cloud  with  charge  density  —qp(r),  where  p(r)  =  Yu=i  l!/h(r)l2-  The  interaction 
potential  between  an  electron  at  r  and  all  other  electrons  is  then  approximated  as  an  integral  over 
P(r)  (77, 12). 


-v  y  zl 

^nE0  Z_i  |r  —  r, 

j= i 1  j 


P(r') 

|r-r'| 


dV  -  Fe_e(r)  . 


(12) 


The  integration  is  over  all  possible  values  of  r',  that  is,  over  all  space. 


With  this  approximation,  the  one-electron  Schrodinger  equation  may  be  solved  iteratively,  as 
illustrated  by  the  pseudocode  in  figure  6.  One  begins  with  an  initial  guess  for  xpf r),  calculates 
p(r)  and  then  Ve_e  and  Text, eft  from  that  guess,  solves  the  one-electron  Schrodinger  equation 
using  the  estimated  Fext,eff?  and  then  checks  if  the  value  of  ipi(r)  and  the  initial  guess 
approximately  match,  given  a  certain  tolerance.  If  not,  the  more  recently  calculated  estimate  of 
xpi(r)  becomes  the  starting  guess  for  the  next  iteration.  (This  iterative  process  could  be  refined 
to  make  it  more  numerically  stable  by,  for  example,  using  a  mixture  of  previous  guesses  as  the 
guess  for  the  next  iteration  [27].)  This  is  called  a  self-consistent  approach. 

The  pseudocode  shown  in  figure  6  does  not  specify  the  means  of  solving  the  eigenproblem  in 
each  iteration.  It  is  possible,  for  example,  to  discretize  the  problem  via  a  finite-difference 
scheme  (22)  or  by  finite  elements  (23).  More  commonly,  though,  the  eigenproblem  is  solved 
through  the  means  shown  in  figure  7,  where  the  one-electron  wavefunction  is  expanded  in  a 
linear  combination  of  basis  functions,  such  as  plane  waves,  atomic  orbitals,  or  Gaussians  (27), 
and  a  matrix  equation  is  obtained  by  substituting  the  expansion  into  the  eigenproblem, 
multiplying  through  by  the  complex  conjugate  of  one  of  the  basis  functions  and  integrating.  The 
resulting  matrix  equation  is  a  generalized  eigenvalue  problem  that  may  be  solved  numerically. 
When  Fext.eff  depends  on  the  one-electron  wavefunctions,  the  latter’s  basis  function  expansions 
are  substituted  into  it  as  well. 
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for  /  =  1  to  n  do 

Make  initial  guess  y/f(r)  for  wavefunction  y/,(r). 
end  for 
Let  m  =  0. 

Set  tolerance  to  tol. 

repeat 

Increment  m  by  1 . 

Let  p"'(r)  =  (r)|2. 

Let^Vcrr(r)  =  K-x,(r)  +  l/'"e(r). 

for  i  =  1  to  n  do 


S°lve  (v2^+V"(r)J  VC(r)=£/V^"(r). 

end  for 


until  err(m,/w  —  1)  <  tol 

where  err(m,m-  1)  is  |y/f'(r)  -  (**)|,  |pm(r)  -pm_l(r)|, 

or  |V"'(r)  — Vm-I(r)|  for  all  r,  i 


Assume  y/,(r)  «=  yf(r),  p(r)  «  p"'( r),  etc. 


Figure  6.  Example  pseudocode  for  iteratively  solving  the  one-electron 
Schrodinger  equation  using  the  Hartree  approximation. 
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Figure  7.  Pseudocode  for  transforming  a  general  one-electron 

Schrodinger  equation  into  a  matrix  equation  that  can  be 
solved  numerically.  The  superscript  “f”  indicates  the 
complex  conjugate. 

The  approximation  in  equation  12  is  fairly  drastic  because  it  implies  that 
lF0({r}n)  =  nf=i<Mri),  which  violates  the  antisymmetry  condition  for  the  many-electron 
wavefunction  in  equation  3.  The  Hartree-Fock  method  (24)  corrects  this  particular  problem  by 
assuming  that  T0({r}n)  is  a  Slater  determinant  of  the  one-electron  wavefunctions. 

tf'ife)  •••  n) 

^o({r}J  =  ^2(ri)  n) 

^n(r r)  ^n(r2)  -  i/>n(rn) 

An  exchange  of  two  coordinates  implies  a  switching  of  two  columns  of  the  determinant,  which 
causes  a  change  in  the  sign  of  vF0({r}n).  As  a  consequence  of  this  form  of  the  many-electron 
wavefunction,  e(r)  becomes  (12) 


where  the  superscript  “f”  indicates  the  complex  conjugate  and  8S.S.  equals  one  if  sL  =  Sj,  and 
zero  otherwise  where  and  s;  are  the  spins  of  electrons  i  and  j.  Pseudocode  for  the  Hartree- 
Fock  approximation  is  shown  in  figure  3.  One  may  solve  the  resulting  eigenproblems  from  this 
algorithm  through  the  method  shown  in  figure  8.  While  in  the  original  Hartree  approximation, 
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for  i  =  1  to  n  do 

Make  initial  guess  V/f(r)  •or  wavefunction  y/,(r). 
end  for 
Let  m  =  0. 

Set  tolerance  to  tol. 


repeat 

Increment  m  by  1 . 


Le>Pra(r)  =  IIV/r'(r)|2. 


i=  1 


Let  V”e(r)  =  -f-  [ f^dV. 
4 KEoJ  r-r' 


for  i  =  1  to  n  do 


Add-- 


a2  n  /•  [t/.',-|(r)lV"“  (r') 

—  tf[j  |r_^  toVc"ic(r)^'(r) 


47t£o 


Let  i>"Vcfr(r)K'(r)  =  [Vex,(r)  +  ^_e(r)]vC(r). 

SolVe  (^  +  ^'".-eff(r))  VO)  =  E;^"( r). 

end  for 


until  err(m,  m  —  1 )  <  tol 

where  err(m,m-  1)  is  | V/T'(r)  —  VT*- 1  (r)U  |p'"(r)  —  p'"—  1  (r)|, 
or  |Vm(r)  — V"”_l(r)|  for  all  r,  / 

Assume  i//,-(r)  ss  !//"'( r),  p(r)  «  p'"( r),  etc. 


Figure  8.  Example  pseudocode  for  iteratively  solving  the  one-electron 
Schrodinger  equation  using  the  Hartree-Fock  approximation. 

Ve_e  is  simply  a  function  of  r  that  is  multiplied  by  t/^( r),  in  the  Hartree-Fock  approximation, 
Ve_e  becomes  an  integral  operator,  and  one  that  involves  n  integrations  for  each  one-electron 
wavefunction  xpi( r).  Taking  the  Pauli  principle  into  account  thus  entails  increased 
computational  expense. 

The  Kohn-Sham  equations  of  density  functional  theory  (DFT)  (25)  are  different  from  the 
previously  mentioned,  effective  one-electron  Schrodinger  equations  because  they  are  essentially 
a  mathematical  tool  to  find  the  total  energy  and  electron  density  at  the  ground  state.  In  principle, 
the  total  energy  of  the  system  is  a  functional  of  the  electron  density  E[p\.  The  value  of  p(r)  that 
minimizes  E  [p]  is  the  ground-state  electron  density,  and  the  minimum  value  of  E  [p]  is  the 
ground-state  total  energy  Ejot  0 .  However,  E  [p]  cannot,  in  general,  be  expressed  as  an  explicit, 
analytic  functional  of  the  density  that  can  be  directly  minimized  (27).  To  work  around  this, 
Kohn  and  Sham  ( 26)  posited  a  fictitious  auxiliary  system  of  noninteracting  electrons  whose 
ground-state  charge  density  and  total  energy  are  identical  to  that  of  the  real  n-electron  system. 
For  this  auxiliary  system,  p(r)  =  £f=i|i/h(r)|2,  and 
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(15) 


i4-e(r)  = 


47T60 


J l~^d3r'  +  Vxc[p] , 


where  Vxc  [p]  is  a  functional  derivative  with  respect  to  p  of  the  exchange-correlation  energy 
Exc  [p],  which  takes  into  account  the  part  of  the  electron-electron  interactions  not  taken  into 
account  by  the  first  term  of  t4-e(r)-  The  Kohn-Sham  equations,  then,  are  solved  iteratively  with 
the  self-consistent  approach  used  to  solve  the  one-electron  Schrodinger  equations  in  the  Hartree 
and  Hartree-Fock  approximations.  Once  the  equations  are  solved,  the  ground-state  total  energy 
may  be  determined  as  follows  (25,  26): 


pn  _  yn  p 
cTot,0  —  Zu=l 


4ne0 


II 


p(r)p(r') 


d3r  d3r'  —  / p(r)Vxcdfr  +  Exc[p]  . 


(16) 


In  principle,  the  Kohn-Sham  scheme  can  yield  the  exact  values  of  the  ground-state  density  and 
total  energy,  though  in  practice  Exc  is  not  known  exactly  and  must  be  approximated.  However, 
because  the  Kohn-Sham  equations  are  one-electron  Schrodinger  equations  for  a  Active,  auxiliary 
system  of  noninteracting  electrons,  there  is  no  rigorous  justification  to  treat  the  eigenvalues  £) 
that  come  from  solving  these  equations  as  addition  or  removal  energies  of  electrons  of  the  true 
n-electron  system.  This  means,  for  example,  that  Ek+1  —  Ek  may  not  even  be  a  good 
approximation  for  the  energy  needed  for  an  electron  to  be  excited  from  energy  level  k  to  k  +  1. 
Often,  the  band  structure  estimated  from  treating  Kohn-Sham  eigenvalues  as  if  they  were  one- 
electron  energies  will  capture  qualitative  trends  but  underestimate  band  gaps  (13). 

There  are  several  other  ab  initio  methods,  which  will  be  mentioned  only  briefly  here  since  they 
are  largely  outside  the  scope  of  this  survey.  The  configuration  interaction  method  builds  upon 
the  Hartree-Fock  method  and  expands  the  many-electron  wavefunction  into  a  linear  combination 
of  Slater  determinants  rather  than  a  single  one  as  in  equation  1  (11).  Both  GW  and  Bethe- 
Salpeter  are  variations  of  many-body  perturbation  theories.  The  GW  method  is  a  generalization 
of  the  Hartree-Fock  approximation  that  accounts  for  dynamic  Coulomb  screening  and  can 
represent  charged  excitations  to  yield  band  structures  that  are  a  closer  match  to  experiment  (13). 
The  Bethe-Salpeter  equation  is  derived  from  the  so-called  four-point  Dyson  equation,  which 
therefore  can  represent  excitonic  effects  and  may  be  suited  for  accurate  estimations  of 
photoabsorption  spectra  (13, 15).  And,  unlike  the  original  DFT  for  time-independent  systems, 
time-dependent  DFT  is  capable  of  yielding  excitation  energies  and  photoabsorption  spectra  (15, 
27).  However,  like  regular  DFT,  it  is  a  reformulation  of  the  many-body  quantum  mechanical 
equation  in  terms  of  electron  densities  and  is  dependent  on  an  available  approximation  for  the 
exchange-correlation  energy  Exc ■  This  can  lead,  for  example,  to  difficulties  in  detennining  the 
excited  states  of  nonmetallic  bulk  solids  (27). 

The  examples  of  ab  initio  methods  briefly  outlined  have  several  points  in  common.  First,  they 
are  free  of  empirical  parameters  aside  from  fundamental  constants  such  as  h  or  me.  Second,  the 
operator  that  each  uses  to  account  for  interactions  between  electrons,  (4-e(r)> 
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depends  on  the  very  one -particle  wavefunctions  xpi  (r)  one  is  seeking  to  solve.  In  principle,  one 
can  solve  the  effective  one -particle  Schrodinger  equation  self-consistently,  that  is,  start  from  trial 
values  of  i/h(r)  and  iterate.  In  practice,  the  Hartree-Fock  and  Kohn-Sham  equations  are  usually 
solved  this  way.  Third,  any  effects  of  strain  can  be  taken  into  account  entirely  through  changes 
in  the  positions  of  the  atomic  nuclei  (R}M,  which  manifest  through  changes  in  the  external 
potential  l/ext.  Finally,  these  methods  are  all  very  computationally  expensive  and  are  suitable 
primarily  for  very  small  systems.  For  example,  as  a  rule  of  thumb,  given  the  current  state  of  the 
art,  DFT  is  usually  feasible  for  systems  of  up  to  a  few  hundred  atoms  (28),  though  there  are  some 
numerical  algorithms  that  introduce  approximations  that  allow  DFT  to  be  applied  to  systems 
with  thousands  of  atoms  (21).  The  GW  method  is  even  more  expensive  and  generally  suitable 
for  systems  with  no  more  than  a  few  10s  of  atoms  (29,  30).  This  starkly  contrasts  with  the 
scalability  of  the  empirical  methods  discussed  in  the  following  sections. 

2.2  Spin-Orbit  Interaction 

In  much  of  this  survey,  the  effects  of  spin  will  be  ignored  for  the  sake  of  simplicity.  However,  if 
need  be,  to  account  for  spin  in  the  one-electron  picture,  the  wavefunction  may  be  written  as  a 
spinor,  where  (31) 

^ik(r)  =  *l> ik(r)  [J]  +  */4(r)  [J] ,  (IV) 

where  the  symbols  “T”  and  “f”  denote  the  spin-up  and  spin-down  states  of  an  electron.  xpjk( r) 
is,  then,  the  Hermitian  conjugate  of  i/hk(r)-  The  one-electron  Hamiltonian  with  a  spin-orbit 
correction  term  is  (32) 

Hie-  (r,  Woo)  =  U  +  kext,eff(r.  Woo)  +  [VFext,eff  X  a] ,  (18) 


where  c  is  the  speed  of  light  and  a  is  a  vector  operator  whose  elements  are  the  Pauli  spin 
matrices  (31,  32) 


0 

1 


(19) 


The  effect  of  spin-orbit  interaction  on  the  band  structure  of  GaAs  is  shown  in  figure  9.  This 
interaction  introduces  the  split-off  band,  so  called  because  it  splits  away  from  the  valence  band 
maximum.  This  effect  occurs  not  only  in  GaAs  but  is  common  in  semiconductors  (32). 
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(a) 


(b) 


Figure  9.  The  calculated  band  structure  of  GaAs  near  the  T  point  (a)  without  and  (b)  with 
spin-orbit  correction.  The  band  structures  were  estimated  using  the  method  and 
parameters  of  Boykin  et  al.  (JJ). 


3.  Atomistic  Methods 


Empirical  approaches  tend  to  be  less  computationally  expensive  and  thus  scalable  to  relatively 
large  systems,  and  this  holds  true  even  for  the  empirical  methods  that  retain  atomic-level 
resolution,  such  as  the  empirical  pseudopotential  (34)  and  tight-binding  (35)  methods. 
Simulations  with  these  methods  can  be  feasible  for  systems  with  thousands  or  even  millions  of 
atoms.  One  reason  for  their  lack  of  computational  expense  as  compared  to  ab  initio  methods  is 
that  they  do  not  require  iterating  to  self-consistency.  While  these  methods  can  involve 
determining  the  eigenvalues  of  large  matrices,  these  values  only  need  to  be  determined  once  in 
the  course  of  a  calculation.  Another  reason  is  that  the  calculation  of  the  elements  of  these 
matrices  is  comparatively  less  involved  than  those  of  the  matrices  formed  in  the  implementation 
of  the  ab  initio  methods.  For  example,  there  is  no  need  to  integrate  a  functional  of  the  electron 
density  p(r),  which  in  turn  requires  a  sum  over  n  one-electron  wavefunctions.  There  is  a  cost  to 
this  relative  ease  of  computational  expense,  though.  First,  these  methods  require  parameters  that 
are  fit  to  experiment  and/or  previous  ab  initio  results,  and  the  process  of  fitting  these  parameters 
is  not  trivial.  As  an  example,  Boykin  et  al.  (33)  use  a  genetic  algorithm  to  fit  over  30  parameters. 
Second,  these  methods  involve  various  approximations,  such  as  neglecting  certain  terms  or 
assuming  that  various  quantities  may  be  expressed  in  certain  functional  fonns.  Even  the  lack  of 
iterating  until  self-consistency  entails  an  approximation,  one  that  prevents  the  changes  in  the 
positions  of  the  atomic  nuclei  (R)M  from  fully  taking  the  effects  of  strain  into  account,  a  matter 
that  which  will  be  discussed  in  more  detail  in  section  3.3.  Nonetheless,  these  empirical 
approaches  are  useful,  especially  where  ab  initio  methods  are  impractical. 
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3.1  Empirical  Pseudopotential  Method 


The  empirical  pseudopotential  method  (34)  is  an  approach  originally  used  to  determine  the 
electronic  band  structure  of  bulk  crystals  {36,  37)  where  the  effective  potential  is  periodic.  The 
effective  potential  may  be  split  up  as  (38) 

Vext,eff(r,  Woo)  =  VexteffCr,  Woo)  +  VNL  ,  (20) 

where  (R}oo  contains  the  set  of  coordinates  of  all  the  atoms  of  the  infinite  bulk  crystal.  The  tenn 
^e°t,eff(r'  Woo)  is  a  function  of  position  and  is  here  called  the  local  part  of  the  potential,  while 
the  operator  VNL  is  a  correction  term  that  accounts  for  Fext,eff(r<  (Rjoo)  not  being  a  mere  function 
of  position  and  is  called  the  nonlocal  part.  Both  tenns  of  Fext,eff(r>  {Woo)  are  periodic,  and 
^e°t,eff(r'  Woo)  can  be  expanded  as  a  complex  Fourier  series, 

ifteffMRU  =Sc^elf(G)eiCr,  (21) 

where  (10,39) 

n«,„(G)  =i;4dxU(>-.{R}»)e-'Crd3r.  (22) 

Here,  f lc  is  the  volume  of  the  unit  cell  of  the  bulk  crystal  and  the  quantity  G  is  a  reciprocal  lattice 
vector.  Because  the  bulk  crystal  is  periodic,  the  Bloch  theorem  applies,  i.e., 

^tk(r)  —  elk  WkC1")-  Because  of  the  periodicity  of  iqk(r),  it  can  be  expanded  in  a  complex 
Fourier  series,  so  that 

^ik(r)  =  y=I!Gnik(G)  ei(k+G)  r ,  (23) 

where  1/^fJi^  is  a  normalization  pre factor.  Alternatively,  this  may  be  described  as  expanding 
ipiki r)  in  terms  of  plane  waves,  each  with  a  wavevector  k  +  G.  If  one  (1)  substitutes  equation 
23  into  equation  7,  the  one-electron  Schrodinger  equation,  (2)  multiplies  both  sides  of  the 
equation  by  e-1(k+G  )'r  /j£Tc,  (3)  integrates,  and  (4)  substitutes  equation  22,  one  may  obtain  the 
following  matrix  equation  (10): 


I 


'  h2 

2m, 


Ik  +  G|25gg,  +  tteffCG'  -  G)  +  VNL (k,  G',  G) 


^tk(G)  —  ^tk^tkCG') 


where 

f/M(k,G',G)  =  i/nce-|(k+=')''(/Arlei<k+G>''d3r 
and 


(24) 


(25) 
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(26) 


1,  G  =  G' 
.0,  G  A  G' 


Here,  fl*  is  an  arbitrary  volume,  which  for  the  previous  derivation  of  the  matrix  equation  is  fic. 
Given  a  value  of  the  Bloch  wave  vector  k,  equation  24  can  then  be  solved  for  the  eigenvalues 
Eik.  In  principle,  there  are  infinitely  many  values  of  the  reciprocal  lattice  vector  G,  so  to  solve 
the  matrix  equation  numerically,  a  subset  of  the  values  of  this  vector  are  used  that  satisfy  the 
condition 


~~  I  k  +  G  | 2  <  Ecut  ,  (27) 

where  Ecut  is  an  energy  cutoff  (37).  Algorithms  for  finding  reciprocal  lattice  vectors  within  the 
energy  cutoff  are  in  appendix  A. 

If  spin-orbit  coupling  is  taken  into  account,  for  each  wave  vector  k  and  reciprocal  lattice  vector 
G  there  are  two  plane  waves,  one  for  a  spin-up  electron  and  one  for  a  spin-down  electron. 
Equation  24  is  then  modified  to  become 


£GXs[tf(k,G',G)5ss, 


+  Hso  (k,  G',  G,  s',  s )  uik(G,  s)  —  Eikuik( G',  s')  , 


(28) 


where  s,  s'  £  (T,  -1}  , 

H(k,  G',  G)  =  |k  +  G|25gg,  +  FSeffCG'  -  G)  +  VNL(k,  G',  G)  ,  (29) 

AUle 

8ssi  =  S  7  S  .  ,  (30) 

(0,  otherwise 

and  Hso  (k,  G',  G,  s',  s)  is  a  correction  term  for  spin-orbit  coupling.  In  some  formulations,  this 
tenn  is  expressed  analytically  as  a  function  of  k,  G,  and  G’  (38),  while  in  others,  the  spin-orbit 
correction  is  calculated  in  real  space  and  then  numerically  converted  into  Fourier  space  (40). 

Both  the  nonlocal  and  spin-orbit  correction  tenns  may  be  neglected,  as  was  done,  for  example,  in 
the  fonnulation  by  Cohen  and  Bergstresser  (37).  Alternatively,  the  spin-orbit  correction  may  be 
kept  while  neglecting  any  other  nonlocal  corrections,  as  was  done  by  Williamson  and  Zunger 
(41).  Williamson  et  al.  (42)  accounted  for  these  other  nonlocal  corrections  not  through  the  tenn 
VNL,  but  by  multiplying  the  kinetic  energy  operator  p  •  p/2me  (or  equivalently, 
h2  |k  +  G|25GCp/2me,  in  Fourier  space),  by  a  fitting  parameter.  For  the  sake  of  simplicity, 
neither  the  nonlocal  conection  term  nor  the  spin-orbit  tenns  in  the  empirical  pseudopotential 
method  will  be  discussed  further. 
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The  principal  approximation  in  the  empirical  pseudopotential  method  is  to  express 
VexteffCG'  —  G)  as  a  function  of  empirically  determined  parameters.  The  local  part  of  the 
effective  potential  Fe££eff(r<  (R}oo)  is  taken  to  be  a  sum  of  contributions  from  individual  atoms 


(34), 


r£,ff(r.  Moo)  =  zr-i  zJIt""  ”/(<■  -  r,  -  dj) 


(31) 


where  R,  is  the  lattice  vector  of  the  ith  unit  cell  of  the  bulk  crystal,  and  R,  +  d j  is  the  position 
vector  of  the  jth  atom  within  this  cell.  Nper  cell  is  the  number  of  atoms  per  unit  cell.  One  can 
rewrite  the  summation  as 


where 


^per  cell 

^xt,e ff(r,  Woo)  =  ^  Vj( r  -  dj)  , 

7  =  1 


(32) 


Vj(f)  =  ZZiVj(r-Rd.  (33) 

Both  kg°^eff(r,  (R}oo)  and  Vj  (r)  are  periodic  and  can  also  be  expressed  as  complex  Fourier 
series.  The  former  has  already  been  shown  in  equations  21  and  22.  For  the  latter,  one  may  write 

Vj(r)  =ZG^(G)eiGr,  Vj(G)  =±-JaVj(r)e-iGrd3r.  (34) 

Substituting  equation  32  into  equation  22  yields 

Nper  cell 

^ext'effCG)  =  V  [  Vj(r-dj)e~iGrd3r  , 
j ri  iLc 


per  cell 


y  —  f  V:(r)  e*'G*(r+‘*l)d3r 

ft  ncJn/ 

7=1 


and 


Nper  cell 

=  V  e~iGdJVj(G). 

7  =  1 


(35) 
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Vj  (G)  is  an  atomic  form  factor,  and  in  the  empirical  pseudopotential  method,  it  is  either  an 
empirical  parameter  (36,  38,  43)  or  is  expressed  in  terms  of  empirically  determined  parameters 
(10,  34,  37).  Pseudocode  for  an  implementation  of  the  empirical  pseudopotential  method  is 
shown  in  figure  10.  This  implementation  assumes  that  the  one-electron  wavefunctions  have  been 
expanded  in  terms  of  plane  waves,  as  in  equation  23.  For  diamond  and  zincblende  crystals, 
equation  35  is  often  rewritten  in  terms  of  symmetric  and  antisymmetric  form  factors  14(G)  and 
14(G),  rather  than  directly  in  terms  of  the  atomic  form  factors.  For  such  crystals,  Nper  cell  =  2, 

and  if  the  lattice  vector  R  is  taken  to  point  to  the  center  of  the  unit  cell,  then  dj  =  — d  and 
d2  =  d,  where  d  =  (ax  +  a2  +  a3)/8,  with  al5  a2,  and  a3  being  the  lattice  vectors  shown  in 
figure  3c.  Given  Euler’s  formula  el°  —  cos  0  +  i  sin  6,  one  may  rewrite  equation  35  as 

kexteff(G)  -  Pi(G)  +  14(G)]  cos(G  •  d)  +  [[VfG)  -  F2(G)]  sin(G  •  d) 

=  14(G)  cos(G  •  d)  +  il4(G)  sin(G  •  d)  .  1 

The  symmetric  and  antisymmetric  fonn  factors  then  become  fitting  parameters  instead  of  the 
atomic  form  factors  themselves.  The  method  shown  in  figure  10  is  then  modified  by  taking 
14(G)  and  14(G)  as  givens  rather  than  the  atomic  form  factors,  and  by  calculating  Ft!”tceff(G) 
through  equation  36  instead  of  35.  Cohen  and  Bergstresser  (37)  treated  14(G)  and  14(G)  as 
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functions  of  the  magnitude  of  G  in  units  of—,  where  a  is  the  lattice  constant  of  the  crystal,  so 

that  14(G)  =  14  (4“)  and  14(G)  =  14  (■—)•  The  values  of  their  form  factors  for  particular 

values  of  I G  |  were  fit  to  reflectivity  and  photoemission  experiments.  Examples  of  some  of  these 
parameters  for  bulk  materials  Si  and  GaAs  are  shown  in  table  1.  For  Si,  14(G)  =  l/2(G),  which 
means  14(G)  =  0  for  this  material.  If  the  unit  cell  changes  size  or  shape  while  the  relationship 
d  =  (ax  +  a2  +  a3)/8  still  holds,  then 

G  d  =  ^If=1rijbj  •  (ax  +  a2  +  a3)  =  ^=1nt;  and  nt  is  an  integer  .  (37) 


Table  1.  Example  symmetric  and  antisymmetric  form  factors  l4(a|G|/27r)  and  VA(a\G\/2n)  from  Cohen 
and  Bergstresser  (37),  where  a  is  the  lattice  constant  of  the  material  and  the  form  factors  are  in 
Rydbergs. 


Material 

F5(V3) 

14(2) 

14  (Vs) 

*4  (VI) 

14(2) 

14  (Vs) 

Si 

-0.21 

0 

0.04 

0 

0 

0 

GaAs 

-0.23 

0 

0.01 

0.07 

0.05 

0 

The  quantity  G  ■  d,  then,  is  independent  of  strains  that  preserve  the  relative  locations  of  the  atoms 
within  the  unit  cell,  such  as  hydrostatic  strains.  However,  the  form  factors  are  not  independent 
of  strain.  For  example,  Aouina  et  al.  (43)  have  recalculated  the  symmetric  and  antisymmetric 
form  factors  for  AlAs  under  a  range  of  pressures,  as  shown  in  table  2. 
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Atomic  form  factors  V\  (G),  V/G), . . . ,  cc||  (G)  are  given. 

Vectors  di,  cfe,  •  •  • ,  d«  ccll  are  given. 

Wavevectors  k| ,  ki, . . . ,  k^k  are  given. 

£cu,  is  given. 

Nonlocal  potential  part  V/v/.(k.G'.G)  is  assumed  to  be  given. 

Spin-orbit  correction  #so(k,G,,G,/,.s)  is  assumed  to  be  given. 

Boolean  variable  noSpin  is  true  if  there  is  no  spin,  and  false  otherwise. 

for  k  in  k|.k2,...,k,vk  do 

Find  the  Nc,  reciprocal  lattice  vectors  Gj,  G2, . . . ,  that  satisfy 
(fi2/2/we)|k  +  G,|2  <  £cui. 

if  noSpin  then 

Let  matrix  [Hu]  be  of  size  Nq  x  Nq. 
else 

Let  matrix  [Hu]  be  of  size  2Nq  x  2Nq. 

end  if 

for  /  =  1  to  Nc,  do 
G'  =  G, 

for  j  =  1  to  Nc,  do 
G  =  G  j 

'Vjxt  cell 

ir(G'-G)=  £  e_i(G,_G)dmV’m(G'  -  G) 

m=l 

H(k, G'. G)  -  |k  +  G|25gG-  +  Ve^eff(G'  -  G)  +  G', G) 

if  noSpin  then 

Let  /  =  i,J  =  j ,  and  Hu  =  H(k,G',G). 
else 

for  /  in  ],  [  do 
for  s  in  I,  |  do 

Let  7  =  2/  —  1  if  /  =f,  and  let  /  =  2 i  otherwise. 

Let  J  =  2j-\  if  s  =|,  and  let  J  =  2 j  otherwise. 

H,j  =H(k,G',G)8SS'  +  Hso(k,G’,G,s',s) 

end  for 
end  for 
end  if 

end  for 
end  for 

Find  eigenvalues  £,k  and  eigenvectors  ;7,k  of  matrix  [Hu], 

end  for 


Figure  10.  Pseudocode  for  empirical  pseudopotential  method  where  the  one- 
electron  wavefunctions  are  expanded  in  a  series  of  plane  waves. 
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Table  2.  Symmetric  and  antisymmetric  form  factors  Vs(a\G\/2n)  and  Va(cl\G\/2tz)  for  AlAs  from  Aouina 
et  al.  (43),  calculated  for  several  pressures.  Pressure  is  in  kilobars.  Lattice  constant  a  is  in 
angstroms  and  the  form  factors  are  in  Rydberg. 


Pressure 

0 

30 

60 

90 

120 

Fs(V3) 

-0.212694 

-0.212931 

-0.212725 

-0.211967 

-0.210679 

PS(V8) 

0 

0 

0 

0 

0 

Fs(VTT) 

0.09275 

0.099816 

0.106973 

0.114209 

0.121637 

Va(V3) 

0.068833 

0.068074 

0.066888 

0.065781 

0.064803 

VA(  2) 

0.05 

0.05 

0.05 

0.05 

0.05 

vA(y/TT) 

-0.0075 

-0.0075 

-0.0075 

-0.0075 

-0.0075 

Lattice  constant  a 

5.6611 

5.6012 

5.5406 

5.4875 

5.4403 

Caruthers  and  Lin-Chung  (44)  work  with  the  atomic  form  factors  directly  in  their  work  on  a 
GaAs-AlAs  superlattice  consisting  of  monoatomically  thin  layers  of  Ga,  Al,  and  As.  The  unit 
cell  of  this  superlattice  is  shown  in  figure  1 1  and  consists  of  four  atoms.  For  this  case,  equation 
35  evaluates  to 

kSeffCGO  =  VGa(G)e~iGd^  +  FA1(G)e-iGdAi  +  vAs(G)[e~iGd^  +  (38) 

where  dGa,  dA],  dAsl,  and  dAs2  are  the  vectors  indicating  the  positions  of  the  Ga,  Al,  and  two  As 
atoms  within  the  unit  cell.  Caruthers  and  Lin-Chung  take  dGa  to  be  zero.  Like  Cohen  and 
Bergstresser,  they  take  the  atomic  form  factors  to  be  functions  of  a  |  G  |  / 2zr  and  treat  them  as 
empirical  parameters.  A  sample  of  these  is  shown  in  table  3.  The  values  of  these  form  factors 
are  fitted  in  order  to  reproduce  the  electronic  band  structures  of  GaAs  and  AlAs.  These  form 
factors  are  no  less  strain-dependent  than  those  of  Cohen  and  Bergstresser  or  Aouina  et  al. 
However,  in  principle,  there  is  nothing  in  the  scheme  of  Caruthers  and  Lin-Chung  that  limits  dA1, 
dAsi,  and  dAs2  to  their  equilibrium  values,  and  they  may  change  in  response  due  to  distortions 
from  strain.  When  applying  the  empirical  pseudopotential  to  more  general  nanostructures,  the 
atoms  of  the  structure  may  be  regarded  as  being  within  a  large  unit  cell,  with  d1;  d2,  d3, ...  as  the 
positions  of  the  atoms  forming  the  nanostructure. 
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Figure  11.  Unit  cell  of  GaAs-AlAs 
superlattice  from 
Caruthers  and  Lin-Chung 
(43).  a  is  the  lattice 
constant  of  GaAs.  The 
figure  was  created  with 
Jmol  (45). 

Table  3.  Example  empirical  parameters  from  Caruthers  and  Lin-Chung  (44),  where  a  is  the 
lattice  constant  of  the  material,  and  the  atomic  form  factors  are  in  Rydberg. 


a\G\/2n 

TGa(a|G|/27r) 

TA1(a|G|/27r) 

TAs(a|G|/27r) 

0 

-0.1114 

-0.1056 

-0.1619 

1 

-0.0750 

-0.0930 

-0.1250 

V2 

-0.0560 

-0.0760 

-0.0970 

Rather  than  use  the  values  of  atomic  form  factors  at  particular  values  of  a\G\/2n  as  fitting 
parameters,  one  may  instead  express  the  form  factors  as  continuous  functions  that  contain 
empirical  parameters.  For  example,  Wang  and  Zunger  ( 46)  take  the  form  factors  to  be 


b(G)  =  and  G  =  |G| , 

a]3e  l4  -1 


(39) 


where  a;1  through  are  fitting  parameters.  Similarly,  Mader  and  Zunger  ( 47)  express  the  form 
factors  with  the  following  function, 


Vj(G)  =  Qj[l  +  fje  W^UcLije  Ci^G  bi^  ,and  G  —  |G| , 


(40) 


where  the  fitting  parameters  are  Slj,  fj,  (3j,  a i;-,  and  cLj.  They  also  fit  these  parameters  for  a 

range  of  unit  cell  volumes,  so  that  the  fitting  parameters  themselves  have  the  desired  property  of 
being  independent  of  at  least  the  hydrostatic  strain.  Here,  the  reciprocal  lattice  vector  G  is  not 
normalized  by  a  factor  of  a/2n,  so  changes  in  a  due  to  strain  are  reflected  in  G.  Also,  since 
these  form  factors  are  continuous  functions  of  G,  the  allowed  values  of  G  are,  at  least  in 
principle,  arbitrary. 
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The  large  cell  containing  a  nanostructure  (or  a  representative  unit  cell  of  it)  may  be  regarded  as  a 
supercell  composed  of  smaller  cells  at  least  approximately  resembling  those  of  a  bulk  crystal. 
^ext,eff(r'  (R)oo),  then,  may  be  written  as 


kSefffr.lRU 


Y1  co  yWcell  vWPer  celltt) 
2ii=l2jy=i  Lk=1 


Vj,K( r  -  R i-  R;  -  dK(/))  » 


(41) 


where  R,  is  the  lattice  vector  pointing  to  the  ith  supercell,  Rj  +  R;  is  the  vector  pointing  to  small 
cell  j  within  supercell  i,  and  d K(j)  is  the  relative  position  of  atom  k  within  small  cell  j.  As  long 
as  the  wavefunction  is  expanded  in  terms  of  periodic  functions — even  if  those  functions  are  not 
plane  waves — a  summation  over  all  the  supercells  tiled  out  to  infinity  is  needed.  The  unit  cells 
are  not  all  alike,  and  may  have  different  kinds  of  atoms  within  them.  This  is  why  Nper  celi, 
d K(j),  and  VjK  are  denoted  here  as  depending  on  the  index  for  each  small  cell  j.  These  cells  may 
even  be  distorted  due  to  strain.  The  preceding  summation  may  also  be  rewritten  as 

V^effCr.  Woo)  =  zr=i iJLT  Xa  Wa(j)va(r  —  Rj  —  R;  -  da)  ,  (42) 

where  a  is  one  of  the  atom  types  in  the  system,  e.g.,  Ga,  In,  As,  etc.,  and  da  is  now  the  relative 
position  of  an  atom  of  type  a  within  one  of  the  small  cells  of  the  system.  Since  a  given  atom 
type  may  not  be  in  a  particular  unit  cell,  the  summation  includes  a  weighting  factor  Wa  (j )  that  is 
0  if  atom  a  is  not  within  small  cell  j  (48,  49).  If  there  is  no  strain  in  the  system,  then  it  is  1  if 
there  is  an  atom  at  da  (49).  Substituting  equation  33  into  the  previous  summation  yields  the 
decomposition  of  Fext,eff(r>  (R}oo)  found  in  Wang  et  al.  (48)  and  Wang  and  Zunger  (49): 

Nce\\ 

^xt,eff(r,  Woo)  =  X  Z  W«WV«(r  ~  -  d«)  ■  (43) 

j= 1  a 

In  the  formulations  by  Kim  et  al.  (50)  and  Williamson  et  al.  (42),  the  weighting  factor  from  the 
previous  equations,  Wa(j),  becomes  strain-dependent  and  has  the  assumed  form 

Wa(j)  =  l  +  yaTr[e(j)],  (44) 

where  ya  is  a  fitting  parameter  and  Tr[e(/)]  is  the  trace  of  the  strain  tensor  at  cell  j.  Williamson 
and  Zunger  (41)  have  found  that  without  an  explicitly  strain-dependent  pseudopotential,  changes 
in  the  band  gap  Eg  due  to  strain  may  be  fit  to  experiment,  but  not  the  changes  in  the  conduction 
band  minima  and  valence  band  maxima  themselves. 


While  the  empirical  pseudopotential  method  has  traditionally  used  a  Fourier  expansion  of  the 
one-electron  wavefunction,  other  basis  set  expansions  can  be  used.  For  example,  in  the  linear 
combination  of  bulk  bands  (LCBB)  form  of  the  method  (48,  49), 


'Pit r)  =  laEmZk^k^mkCr).  and  xp^k( r)  =  eikr<k(r)  ,  (45) 
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where  Cfk(J  is  an  expansion  coefficient,  i/'mkC1')  's  the  one-electron  wavefunction  with  band 
index  m  for  a  bulk  material  o  (e.g.,  GaAs  or  InAs),  and  the  values  of  k  are  a  sampling  of  points 
within  the  Brillouin  zone  of  the  bulk  material  o.  Even  though  the  one-electron  wavefunction  is 
no  longer  directly  expanded  in  a  basis  of  plane  waves,  the  previously  shown  atomic  form  factors 
expressed  in  Fourier  space,  such  as  those  in  table  3  or  equations  39  and  40,  may  still  be  used. 

This  is  because  it^k(r)  is  still  expanded  in  a  plane-wave  basis,  so  that 

«kW  =  7=S  Gfi"k(G)e'<k+G>'\  (46) 

where  fi"  is  the  volume  of  the  primitive  cell  of  bulk  material  o.  The  expansion  coefficients 
u"k(G)  may  be  detennined  through  applying  the  traditional  plane -wave  based  empirical 
pseudopotential  method  to  a  unit  cell  of  material  er  for  a  chosen  value  of  Ec ut,  following  the 
method  illustrated  in  figure  10.  For  the  case  where  the  nonlocal  part  of  the  potential  is  zero  and 
there  is  no  spin-orbit  correction,  the  matrix  equation  to  be  solved  is 


where 


V  V  V  J-f  ®  ® 

L.O  Aim  Aik  nm'k'mk^mka 


FT1  ■  ■ 
ciLm  k< 


(47) 


Tjo  a 
nm'k'mk 


h2a 

2mP 


^  “mV(G,)“mk(G)lk  +  G|2<W<W 

G,G' 


+  ^“mk'(G^mk(G) 


wcell 

x  V  Va(k'  +  G'  -  k  -  G)e-id"“  (k+G-k'-G')  V  VKa0')ei(k_k')  Rf  (48) 

a  j=l 

lb  is  the  volume  of  the  simulation  domain,  Va  is  the  atomic  form  factor  expressed  in  Fourier 
space,  and  Wa(j)  is  the  weighting  factor  of  equation  43.  There  is  also  a  variant  of  the  FCBB 
pseudopotential  method  called  the  strained  linear  combination  of  bulk  bands  (SFCBB),  where 
the  basis  function  i/>^k( r)  is  A//(r)elk'ru£lk(x-1(r)),  where  x  1  is  an  inverse  deformation  map 
that  maps  coordinate  r  in  the  deformed  system  onto  its  counterpart  r0  in  the  undeformed  system, 
and  d3r°  =  /  (r)d3r.  In  a  sense,  the  FCBB  and  SFCBB  methods  can  be  described  as  indirect 
plane  wave  expansions.  A  one-electron  wavefunction  of  the  whole  system,  i/h(r)>  is  expanded 
in  terms  of  the  one-electron  wavefunction  of  its  constituent  bulk  materials,  (V;mk(r)}?  which  are 
then  expanded  into  plane  waves.  This  indirect  expansion,  though,  means  that  the  dimensions  of 
the  matrix  [ff^/k'mk]  are  independent  of  the  number  of  plane  waves  used  to  expand  ^mkM- 

The  overall  speed  of  the  empirical  pseudopotential  method  can  be  illustrated  from  a  couple  of 
examples.  A  parallel  code  implementing  the  traditional  plane-wave  version  of  the  method  ran  on 
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a  Cray-T3E900  cluster.  In  this  implementation,  the  one-electron  wavefunction  was  taken  to  have 
the  plane  wave  expansion  shown  in  equation  23.  The  atomic  form  factors  were  from  Wang  et  al. 
(51),  and  the  energy  cutoff  was  set  to  Ecut  =  5  Ry.  The  code  was  used  to  simulate  a  system 
composed  of  an  InAs  quantum  dot  embedded  in  a  GaAs  matrix.  For  a  system  containing 
250,000  atoms,  this  code  took  about  20  h  to  run  on  128  processors  of  the  cluster  (52,  51).  A 
serial  code  implementing  the  SLCBB  version  of  the  empirical  pseudopotential  method  (53)  was 
able  to  simulate  the  same  250,000-atom  system  in  about  10  h  on  an  IBM  595  workstation,  and 
the  eigenvalues  determined  by  this  code  were  within  0.8%  of  those  calculated  by  the  parallel 
plane-wave  code.  More  details  of  the  SLCBB  simulation,  such  as  the  sampling  of  the  Brillouin 
zone,  the  bulk  bands  chosen,  etc.,  are  in  Wang  and  Zunger  (49).  Both  codes  used  the  same 
method  to  detennine  eigenvalues  of  a  matrix,  the  folded-spectrum  method  (46),  which  scales 
linearly  with  the  dimension  of  the  matrix.  The  scalability  of  the  empirical  pseudopotential 
method  is  also  partially  determined  by  the  choice  of  basis  functions.  Now,  both  the  plane- wave 
and  LCBB/SLCBB  forms  of  the  empirical  pseudopotential  method  require  a  sum  over  the  atoms 
in  the  simulation  for  each  matrix  element,  at  least  if  spin  is  not  taken  into  account.  In  the 
pseudocode  shown  in  figure  4  for  the  plane-wave  method,  this  is  apparent  in  the  calculation  of 
'/cxt,cff(r'  Woo),  where  there  is  a  summation  over  all  the  atoms  in  the  unit  cell,  and  it  is  implied 
in  equation  48  for  the  LCBB  method,  where  there  is  a  summation  over  the  jVcen  small  cells  that 
compose  the  supercell  of  the  system  and  a  summation  over  the  atom  types  a  in  the  system. 


However,  in  the  traditional  plane-wave  empirical  pseudopotential  method,  another  scalability 
concern  is  that  the  size  of  the  matrix  ]  whose  eigenvalues  are  to  be  found  is  determined  by 
the  number  of  plane  waves  NG  ,  whose  reciprocal  lattice  vectors  satisfy  an  energy  cutoff  criterion 
(i.e.,  equation  27).  This  number  may  be  roughly  estimated  from  equations  A-l  and  A-2  in 
appendix  A  to  be 

-  W“)(n S“>  W)  =  (M 


This,  in  turn,  is  proportional  to  the  volume  H  of  the  unit  cell  of  the  simulated  system.  Since 
•  a j  —  2nSij,  the  magnitude  of  a  primitive  reciprocal  lattice  vector  is  |bj  =  2n/ (laj  cos  L 
where  6^  is  the  angle  between  b,  and  its  corresponding  primitive  lattice  vector  at.  Therefore, 


Nr. 


fr T3  lail  cos 9{\  I 
\[li  =  1  2n  )\ 


2meEcut 


h2 


oc  It 


(50) 


because  It  oc  nLilat  I  •  This  scaling  with  the  volume  holds  regardless  of  whether  the  volume  is 
filled  with  atoms,  so  it  would  apply  if  the  plane-wave  pseudopotential  method  were  applied  to, 
for  example,  a  cluster  of  atoms  at  the  center  of  an  otherwise  empty  unit  cell  that  was  made  large 
enough  to  prevent  the  cluster  from  interacting  with  its  periodic  images.  In  the  LCBB  and  SLCBB 
forms  of  the  empirical  pseudopotential  method,  though,  the  number  of  plane  waves  used  to 
expand  each  basis  function  xpmk( r)  is  independent  of  the  volume  or  even  the  number  of  atoms  in 
the  simulation.  The  LCBB  and  SLCBB  methods  are  clearly  more  scalable  than  the  plane -wave 
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version  of  the  empirical  pseudopotential  method,  but  the  matrix  elements  in  the  former  two 
methods  are  more  computationally  expensive  to  calculate,  so  such  methods  may  not  be  ideal  for 
smaller  systems.  Also,  application  of  the  LCBB  and  SLCBB  methods  requires  a  careful  choice 
of  the  bulk  band  indices  and  Brillouin  zone  sampling  (i.e.,  m  and  k  in  equations  45-48).  In 
short,  there  are  tradeoffs  between  these  two  forms  of  the  empirical  pseudopotential  method. 

3.2  Slater-Koster  Tight-Binding  Method 

In  tight-binding  methods  {54-59),  the  one-electron  wavefunction  is  expanded  in  terms  of  either 
atomic  orbitals  or  functions  with  the  same  symmetries  as  atomic  orbitals.  For  a  finite  system  of 
M  atoms,  one  may  write  (56,  57) 


i/h(r)  =  Xm=i£oCo^0ozm(r-Rm)  ’ 


(51) 


where  o  is  a  type  of  orbital  symmetry,  i.e.,  o  G  {s,  px,  py,  pz,  dxy, ... },  C^m  is  an  expansion 
coefficient,  and  4>ozm  (r  —  Rm)  is  the  atomic  orbital  or  orbital-like  function  centered  at  position 
vector  Rm,  with  symmetry  o  for  an  atom  whose  species  is  indicated  by  atomic  number  Zm. 

4>ozm(r )  is  real,  so  <plzm (r)  =  <fi0zm (r) •  Also,  <p0zm (r)  is  localized  and  therefore  decays  to  zero 
as  |r|  — >  oo.  For  an  infinite  periodic  system  that  satisfies  the  Bloch  theorem,  the  wavefunction  is 
expanded  instead  as  {54,  57,  60) 

^per  cell 

^Co(m0okm(r)  ,  (52) 

m=  1  o 

where  Nper  cell  is  the  number  of  atoms  per  unit  cell,  is  an  expansion  coefficient,  and 

(pokm  (r)  is  a  Bloch  sum,  defined  such  that 

0okm (r)  =  lim^^oo  eik'(R;+dm)0oZm(r  -  R;  -  dj  .  (53) 


0ik(r)  =  ^ 


Here,  Acell  is  the  number  of  unit  cells  in  the  system,  Ry  is  the  lattice  vector  pointing  to  unit  cell  j, 
and  R;  +  dm  is  the  position  of  atom  m  in  cell  j.  If  one  substitutes  the  tight-binding  expansion 

for  a  finite  system,  i.e.,  equation  51,  into  the  one-electron  Schrodinger  equation  7,  multiplies 
both  sides  of  the  equation  by  <p0'Z  ,  (r  —  R™').  and  integrates  over  all  space,  then  one  obtains 

m 

the  following  matrix  equation: 
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m=l  o 
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o'm'  om 


{-‘om 


M 


m=  1  o 


(54) 


where 
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(55) 


Ho'm'om  =  I  4>o'z  ,  ~  Rm0  ^le“  0oZm  (r  “  Rm)d3r 

m  "L 

and 

So'm'om  =  J  0o'zm,(r  -  Rm0  0ozm(r  -  Rm)d3r .  (56) 

Matrix  equation  54  is  a  generalized  eigenvalue  problem  that  may  be  solved  for  EL.  S0'm'0m  is 
called  the  overlap  matrix,  and  if  it  were  diagonal,  that  is,  if  S0>m>0rn  —  8m'm80'0,  then  the 
generalized  eigenvalue  equation  would  become  an  ordinary  matrix  eigenvalue  equation,  which  is 
less  computationally  expensive  to  solve  than  the  generalized  equation.  However,  if  4>ozm  (r 
— Rm)  is  a  true  atomic  orbital,  that  is,  a  one-electron  wavefunction  with  symmetry  o  of  an  atom 
with  atomic  number  Zm,  then,  in  general,  S0'm'0m  A  8m'm80i0.  Since  orbitals  centered  by  the 
same  atom  are  orthogonal  to  one  another,  S0imom  =  80<0  still  holds.  Also,  if  atoms  m!  and  m 
are  sufficiently  far  apart,  then  S0'm'0m  »  0  because  the  orbitals  are  localized.  If  (p0zm  (r  ~  Rm) 
is  a  Lowdin  orbital,  that  is,  given  a  set  of  true  atomic  orbitals  {4>0zm  (r  —  Rm)},  0ozm(r-Rm) 
is  such  that 

0„z„(r  -  Rm)  =  -  M . 

o'  m 

~  0oZm(^"  Rm) 


'  9  ^  y.goWom  ^o  o^m'm)<^o  Zm’(r  ^m  ) 


(57) 


and 


8 o'm'om  f*Po'z  t  (T  RmO  0oZm(^"  Rm)d  If 


(58) 


where  a  superscript  of  —  1/2  indicates  the  inverse  square  root  of  a  matrix,  then  it  can  be  shown 
that  S0'm'0rn  becomes  diagonal  and  equals  80f08mfrn  (61). 


The  tight-binding  matrix  equation  for  an  infinite  periodic  system  is  similar  to  that  for  a  finite 
system,  but  with  Nper  cen  replacing  M,  and 


where 


Nper  cell 

M  /~(ik) 

11  o  m  omkL‘om 

m= 1  o 


Nper  cell 


=  E: 


1  1- 


(tk) 


ik  /  /  L,o'm'omkQm  , 

I 

m= 1  o 


(59) 


^o'm'omk  I  0o'km'  (x)^le  0okm(f)d  T 
and 

^o'm'omk  /  Po'km’  (O0okm(Od  r. 


(60) 

(61) 
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The  integration  here  is  over  a  unit  cell  of  the  periodic  system.  S0im'omk  is  an  overlap  matrix 
similar  to  the  one  for  the  finite  system  S0'm'0m.  In  general,  S0'mornk  =  80/0,  and  if  the  functions 
4>0zm  (r)  in  the  Bloch  sum  0oUm(r)  are  orthogonalized  rather  than  true  atomic  orbitals,  then 

^o'm'omk  ~  ^m'm^o'o’ 

In  the  Slater-Koster  tight-binding  method,  the  Hamiltonian  matrix  elements,  either  H0tmi  om  or 
H o'm'omk'  are  decomposed  into  sums  of  integrals  that  are  neglected,  become  empirical 
parameters,  or  become  linear  combinations  of  special  two-center  integrals  that  are  functions  only 
of  the  distance  between  their  two  centers.  These  special  integrals  are  generally  not  evaluated 
directly  through  integration,  but  rather  are  estimated  through  simple  closed-form  functions  with 
empirical  parameters.  The  decomposition  begins  by  decomposing  the  effective  potential  into  a 
sum  of  per-atom  contributions,  such  as 

M 

hext,eff(r.  Wm)  =  ^  vZk  (r  -  Rfc)  (62) 

k=l 

for  a  finite  system  with  M  atoms,  or 

^percell  oo 

^exteffCr,  Woo)  =  ^  ^  Vzk{y  -  R;  -  d /c)  (63) 

k= 1  j= 1 

for  an  infinite  periodic  system.  (At  this  point,  spin-orbit  coupling  is  ignored.)  If  m  A  m' ,  then 
after  substituting  equations  7  and  62  in  equation  55,  one  may  write  the  Hamiltonian  matrix 
element  as  (62,  63) 


H o'm'  om 


+  \j  <Pozm( r-Rm') 
+  f  0o'zm,  (r-RmO 


p  p 


2  me 

P  P 


+  vz(r-Rm) 


+  vz(r-Rm) 


2me 

vz,  (X  -  Rm0  +  vZm  (r  -  Rm) 


0ozm(r-  Rm)d3r 
0o'z  ,(r-Rm)d3r 


0ozm(r-  Rm)d3r 


+  ^  j  (Po'zmXr-Rm')vzk(X-Rk)fozm(.r-R  m)d3r. 

k±m±m' 

For  the  special  case  where  m  —  m' ,  the  Hamiltonian  matrix  element  is  instead 


(64) 


Hr 


0o 


'Zn,  (r  Rjn) 


p  P 


2  mP 


+  vZm( r-  Rm) 


0ozm(r-Rm)d3r 


+  X  J  &  o' zm(r  ~  Rm)  vZk(r  -  Rk)  <t>0zm  (r  -  Rm)  d3  r . 


k=tm 


(65) 
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At  this  point,  one  may  begin  simplifying.  First,  the  three-center  integrals,  that  is,  the  integrals 
containing  three  factors,  each  centered  around  a  different  atomic  site  (i.e.,  m',  k,  and  m),  are 
treated  as  if  they  were  negligible  (54).  This  removes  the  integrals  in  the  summation  over 
k  A  m  A  m'  in  equation  64.  The  integrals  whose  potential  is  centered  at  site  k,  while  the 
orbitals  are  centered  around  m,  are  neglected  much  as  the  three-center  integrals  are  (64),  which 
removes  the  integrals  in  the  summation  over  k  A  m  in  equation  65.  This  leaves  the  two-center 
integrals,  whose  factors  are  centered  around  m'  or  m,  and  the  onsite  integrals,  all  of  whose 
factors  are  centered  around  m.  As  Slater  and  Koster  themselves  pointed  out,  the  integrals  that 
are  neglected  here  are  not  necessarily  negligible  in  comparison  to  the  two-center  integrals,  but 
they  are  smaller  than  the  two-center  integrals,  and  their  neglect  will  reduce  the  number  of  fitting 
parameters  needed  (54).  Another  simplification  comes  from  noting  that  an  isolated  atom  at  Rm 
with  the  effective  potential  vZm(r  —  Rm)  has  the  following  one-electron  Schrodinger  equation: 


+  vzJy  -  Rm)]  4>oZm  (r  -  Rm)  =  EoZJ)oZJy  -  Rm)  , 


(66) 


where  EoZm  is  the  eigenvalue  corresponding  to  the  orbital  symmetry  o  for  an  atom  with  atomic 
number  Zm.  If  one  denotes  the  two-center  integral  containing  the  average  of  two  atomic 
potentials  as 


■/ 


vz  ,(r  -  Rm/)  +  vz  (r  — Rm) 

0o'z  ,(r  -  Rm0  ,  m  - -0ozm( r  -  Rm)d3r  ,  (67) 


then,  given  the  simplifications,  the  Hamiltonian  matrix  elements  in  equations  64  and  65  become 
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o'm'  om 


E  OZyn+E  oZ  I 

_ m_  r*  |  I/ 

2  ^ o'm'om  '  vo'm'om  • 

Ho' mom  ~  HoZm80'0 


(68) 


The  first  expression  indicates  the  Hamiltonian  matrix  element  for  two  arbitrary  orbitals  o'  and  o 
for  two  different  atoms  m'  and  m,  respectively.  The  second  expresses  the  matrix  element  for  the 
case  where  two  different  orbitals  o’  and  o  are  centered  about  the  same  atom  m.  The  atomic 
eigenvalues  are  fitting  parameters.  Thus,  the  energies  must  be  provided  from  another  source, 
generally  empirical  data  or  through  more  accurate  electronic  structure  ab  initio  approach. 
However,  V0'm'0m  is  not  a  fitting  parameter.  Rather,  provided  that  atoms  m'  and  m  are  not  so 
far  apart  as  to  make  it  negligible,  V0'm'0m  is  decomposed  further  into  sums  of  special  two-center 
integrals  that  depend  only  on  the  distance  Rmm'  between  atomic  sites  m  and  m' .  To  do  this,  a 
change  of  coordinates  r  r  +  Rm’  is  applied  to  V0'm'0m  to  make  explicit  that  it  depends  on 
Rm  R m' ’  th'h  is, 


^o'm'om 


r,  VZm'W+Vzm(r-RmmO  ,  f 

=  1 0o'z  ,(r)  — - ; - 0ozm(r 


R  r 

l^mm 


)d3r . 


(69) 


For  brevity,  the  average  of  the  two  potentials  at  sites  m!  and  m  will  be  denoted  as 
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=  <gm(r)  ■ 


(70) 


VzmXr)  +  Vzm(.r-Rmrn/) 
2 


One  of  the  special  two-center  integrals  centered  around  the  two  nearest  neighboring  sites  m  and 
m!  whose  orbitals  both  have  5  symmetry  is  denoted  here  as  Vs 


r  sm  ,sm, a -> 


V'. 


W 


0s,zm,(r)  ^gm(r)0s,zm(r  -  Rmrn')d3r . 


(71) 


Here,  because  of  the  symmetry  of  the  5-orbitals,  the  direction  of  Rmm'  does  not  matter;  only  its 
magnitude  matters.  Another  special  two-center  integral  involves  the  overlap  of  5-  and  p-orbitals, 
as  shown  in  figure  12.  These  two-center  integrals  for  Rmrn'  along  the  x-,  y-,  and  z-directions  are 

^sm',pm,a  ~  j  0sZm/ (O^avg  (r)0pxZm(r  —  > 

=  |  0szm,(rXgrn(r)0pyzm(r-«mm'ey)d3r  , 

and 

=  |  0szm,  (r)^gm(r)0pzzm (r  -  Rmm'^z) d3r .  (72) 


Figure  12.  Overlapping  between 
s-  and  p-orbitals. 

The  plus  and  minus 
signs  indicate  where 
an  orbital  has  a 
positive  or  negative 
value. 

Other  special  two-center  integrals  involve  pairs  of  overlapping  p-orbitals.  Figure  13  shows  the 
two  types  of  p-orbital  overlap.  Suppose  that  there  are  two  px -orbitals  with  o -type  overlap,  and 
the  vector  connecting  their  centers  points  along  the  x-direction,  that  is,  Rmrn'  —  Rmm'ex,  where 
ex  is  a  unit  vector.  Alternatively,  suppose  that  the  two  orbitals  with  er  overlap  are  of  py- type 
with  Rmm'  —  Rmrnfey,  or  of  pz- type  with  Rmm'  —  where  ey  and  ez  are  unit  vectors 

pointing  along  the  y-  and  z-directions.  By  symmetry,  the  two-center  integral  in  all  of  these  cases 
has  the  same  value  (55): 
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=  /  0pzzm,(rXgrn(O0Pzzm(r  -  ^mm'ez)d3r  . 
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a  overlap  n  overlap 


Figure  13.  Types  of  overlap  between 
p-orbitals.  The  plus  and 
minus  signs  indicate  where 
an  orbital  has  a  positive  or 
negative  value  (55). 


If  the  two  p-orbitals  have  7r-type  overlap  and  Rmm'  =  Rmm'ex,  the  orbitals  must  both  be  either 
py  or  pz,  that  is, 

Vpm' ,pm,n  ~  j  0pyZm/ (O^avg  (O0pyZm(^  —  > 

and 

=  j  0Pzzm,  (r)v£gm(r)c/)pzZm(r  -  Rmm,ex)d3r .  (74) 

Similarly,  if  Rmrn>  —  Rmm'ey,  then  the  orbitals  in  Vprn>  prn  Tt  must  both  be  of  px  or  pz  type,  and  if 
Rmm'  =  ^mm'ez»  then  the  orbitals  in  Vpm'  pm  n  must  both  be  of  px  or  py  type. 


Many  two-center  integrals  turn  out  to  be  zero.  For  example, 

/  0szm,(rXgrn(O0pyZm(r  -  Rmm'ex)d3r)  _  ^ 
/  0szm,(rXgrn(O0pzzm(r  -  Rmm'ex) d3r) 


The  reason  can  be  seen  in  figure  14.  The  interaction  between  the  s -orbital  and  the  positive  lobe 
of  the  p-orbital  is  exactly  canceled  by  a  corresponding  interaction  between  the  s-orbital  and  the 
p-orbital’s  negative  lobe.  For  similar  reasons,  the  two-center  integrals  associated  with  the  kinds 
of  interaction  between  p-orbitals  shown  in  figure  15  are  also  zero  ((55),  e.g., 


/  0Pxzm,  (rXgm(r)0pyzm(r  ~  tfmm'e*)d3rj  ^ 
/  0pyzm,(rXgm(r)0pzzm(r  -  Rmm>ex)  d3rj 
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Figure  14.  s-  and  p-orbitals  whose  net  overlap 
is  zero.  The  plus  and  minus  signs 
indicate  where  an  orbital  has  a 
positive  or  negative  value.  The 
overlap  contribution  from  the 
positive  lobe  of  the  p-orbital  is 
exactly  canceled  by  that  in  the 
negative  lobe. 


Figure  15.  Pairs  of  p-orbitals  whose  net  overlap  is  zero.  The 
plus  and  minus  signs  indicate  where  an  orbital  has  a 
positive  or  negative  value.  The  bond  axis  runs 
through  the  points  connected  by  the  vector  Rmm'- 

With  the  special  two-center  integrals  Vsm'  pma,  Vpm> pma,  and  Vprn>  pmn  identified,  the 
decomposition  of  a  general  two-center  integral  V0imi0rn,  where  o',  o  E  {s,  px,  py,  pz],  can  be 
shown  (54,  55).  Let  there  be  a  coordinate  system  x'y'z'  such  that  Rmm'  =  Rmm'e'x,  and  the 
transfonnation  between  from  unprimed  to  primed  coordinates  is  (66) 

x  =  (ex  •  e'x)x'  +  ( ex  ■  e'y)y'  +  (ex  •  e'z)z' , 

y  —  '  &x)x  +  (f!y  ’  "h  (Cy  ’  6z)z  , 

and 

z  =  (ez  •  e'x)x'  +  (ez  •  e'y)y'  +  (ez  ■  e'z)z'. 


(77) 
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The  p-orbitals  can  be  factored  as 


0Pxzm(r)  =«zm(|r|)x, 

0Pyzm( r)  =  RZm(\r\)y, 
and 

0Pzzm(r)  =  Rzm(\r\)z,  (78) 

where  RZm  (|r|)  is  a  function  dependent  on  the  magnitude  of  r  (11).  A  px -orbital,  then,  can  be 


expressed  in  x'y'z'  i 

coordinates  as 

0pxzm(r) 

=  Rzm(.\r\)[(ex  ■ 

6x)x  +  (ex  •  ey)y  +  (ex 

e'z)z'] , 

(79) 

4>pxzjr) 

=  ^Zm(lrl)[(ex 

6x)x  +  (bx  ■  f?y)y  ~b  (&x 

e'z)z'] , 

and 

=  0pizm(r)(ex 

■  ex)  +  (PpyZm(rKex  '  ey)  + 

<j PP'zZm(y)(ex  ■  e'z) , 

(79) 

and  similarly  for  the  py-  and  pz -orbitals.  Accordingly,  a  general  two-center  integral  between  an 
s-orbital  and  a  px -orbital  is 

K sm'pxm  —  j  0sZm»avg  (O0pxZm(r  —  T 

=  /  ■  «l) 

+ 1  0szm,(r)^grn(r)0p;Zm(r  -  iWe*)d3r  (ex  ■  e'y ) 

+ 1  ^s^m'W^g^CO^p^Cr  -  Vei)d3r  (e*  •  e'z ) 

=  7sm'pm,ff(ex  •  e(.)  .  (80) 

A  general  two-center  integral  between  two  p-orbitals  can  be  expressed  as  a  linear  combination  of 
Vpm',pm,a  and  Vpm' ,pm,n •  F°r  example, 

Vpxm'pxm  —  /  0pxZm/ (O^avg  (O0pxZm(^  1", 

and 

—  K  /  (e  /  17 e  •  e;  1  +(g  .e')’l  (81) 

vpm  ,pm,a  ycx  cxJ  ~  v pm  ,pm,n  cy )  '  vcx  czJ  I  •  v  7 

The  directions  of  the  y'  and  z'  axes  are  arbitrary  except  for  the  constraint  that  they  must  be 
orthogonal  to  x'  and  each  other,  and,  accordingly,  the  previous  two-center  integral  can  be 
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rewritten  to  indicate  that  it  does  not  depend  on  e'y  and  e'z  but  rather  on  the  direction  cosines  of 
R  mm',  Ir  =  (e*  •  e'x),  mR  =  (ey  •  e*),  and  nR  —  (ez  •  e*).  Because  of  the  orthogonality  of 
primed  and  unprimed  coordinates,  one  may  write  ( 66) 

(e*  •  e'y)2  +  (e*  ■  e’z )2  =  1  -  (e*  •  e'x )2  =  1  -  l2R  .  (82) 

Expressions  of  general  two-center  integrals  in  tenns  of  lR,  mR,  nR,  and  the  special  two-center 
integrals  have  been  tabulated  in  the  classic  paper  of  Slater  and  Koster  (54).  Implementations  of 
the  tight-binding  method  can  detennine  the  decomposition  of  general  two-center  integrals 
through  a  look-up  table  based  on  this  tabulation. 

The  elements  of  the  overlap  matrix  S0tm'om  in  equation  56  are  similar  in  fonn  to  those  of  the 
general  two-center  integral  V0'm'0m  of  equation  69,  except  that  [vz  ,  (r)  +  vZm  (r  —  Rmm/)]/2 

is  replaced  by  1 .  Therefore,  it  is  possible  to  decompose  S0'm'0m  into  a  linear  combination  of 
integrals  analogous  to  Vsm>sma,  Vsm'pm(T,  Vprn'iPm(T,  Vpm'iPmn,  etc.  That  is,  one  may  define 

(19,  67) 

S sm' ,sm,a  ~  f  $s,Zmi(T)  (ps.ZmtV  T  (83) 


and 


S. sm' ,pm,a  j  4>sZmr(?)<fipxZ7n  C1"  Rmm'^x^6  r  ■> 

=  f  0szm,(r)0pyzm(r-  Rmm'ey)d3r , 
=  /  0szm,(r)0Pyzm(r  -  ^mm'ey)d3r , 


(84) 


pm' , pm, o  anc*  Spm'.pm.n  maY  be  defined  by  replacing  iy^grn(r)  with  1  in  the 
two-center  integrals  Vpmi  pm/J  and  Vpm'  pmjz.  Once  these  integrals  are  defined,  decompositions 


and  so  on.  Integrals  S 
two-center  integrals  V ? 
such  as  the  following  may  be  done: 


and 


S sm'pxm  $  sm'  ,pm,a^x  '  ^x)  sm'  ,pm,a 


Spxm'pxm  $  pm'  ,pm,o  ^ R  $  pm'  ,pm,n 


(1  -  12r)  ■ 


(85) 

(86) 


These  example  decompositions  are  analogous  to  the  decompositions  of  Vsrn'Pxm  and  VPxm'Pxm 

in  equations  80  and  81,  respectively.  As  pointed  out  before,  if  atoms  m  and  m!  are  sufficiently 
far  apart,  then  the  localization  of  the  orbitals  <p0'z  (r)  and  <poZ  (r  —  Rmm ')  entails  that 

m  171 

0.  This  applies  to  the  integrals  Ssm' pma,  Spm> pma,  etc.,  as  well,  and  the  localization 


S o'  m!  om 


of  the  orbitals  also  entails  that  the  two-center  integrals  Vs 


VB. 


sm  ,sm, o  9  r  sm  ,pm 


ff,  etc.  decay  as 


Rmm'  °°- 


Accounting  for  spin-orbit  coupling  means  that  the  number  of  orbital  types  is  doubled,  e.g., 
instead  of  o' ,  o  E  {s,  px,  py, ...  ],  one  has  o' ,  o  £  {s  t,  px  T,  py  T, ,  s  i  px  i,  py  f, ...  ],  where 
again  the  symbols  “T”  and  “f  ”  denote  the  spin-up  and  spin-down  states  of  an  electron.  Let 
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o  =  oS,  where  o  G  {s,  px,  py, ... },  and  S  G  {T,  f}.  An  element  of  a  tight-binding  Hamiltonian 
matrix  with  the  coupling,  then,  may  be  expressed  as  (68) 


Hr 


=  H s 


o  momuS  S 


+  H. 


so 

o'mom  ’ 


(87) 


where  Hdrmdm  is  the  corresponding  Hamiltonian  matrix  element  without  spin-orbit  correction. 

on 

In  the  scheme  of  Chadi  (69),  the  correction  term  H0'mom  is  zero  except  for  the  following  matrix 
elements  (70,  71): 

it SO  _  (  iiSO  ^  _  _rrSO  —  _  (  t/SO  V  =  i  t 

npxlmpylm  ynpylmpxlm )  npx\mpy'\m  ynpytmpxtm )  lAm  > 


HSO  _  (it  SO  H  _  _„SO  _  _(  rrSO  H  _  2 

npzlmpx1m  \rlpx\mpzim)  npz1mpxlm  \rlpximpz/\m)  Am  ■> 


and 


irso  =(mso  V  _  Hso  =(mso  )t  — 

nPylmpz1m  ynpztmpylm )  npy1mpzlm  \npzlmpx1m)  lAm  ■> 


(88) 


where  Am  is  an  empirical  parameter. 

Figure  16  shows  pseudocode  for  the  tight-binding  method  applied  to  finite  system.  It  is  assumed 
in  this  algorithm  that  the  means  of  empirically  estimating  the  special  two-center  integrals 
vsm',sm,a’  Vsm' ,Pm,a>  vPm’.Pm.ai  vVm' ,pm,n>  etc-> is  already  available.  So  that  the  algorithm  shown 
would  fit  on  a  single  page,  an  optimization  was  left  out  that  would  be  used  in  a  more  realistic 
code,  which  would  take  advantage  of  the  fact  that  both  V0im'om  and  S0'm'om  decay  to  zero  for 
large  Rmm'  ■  If  atoms  m'  and  m  are  sufficiently  far  apart  (e.g.,  not  nearest  or,  possibly,  next- 
nearest  neighbors),  then  H o'mom  is  negligible  and  can  be  immediately  set  to  zero  without  an 
explicit  determination  of  V0'm'0m  and  S0imi0m,  which  would  be  approximately  zero. 
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The  orbitals  s,  px,  py,  etc.  for  the  basis  set  are  given.  The  number  of  orbitals  in  the  set  is  Vor b. 
The  number  of  atoms  /Vatoms  is  given. 

Values  of  two-center  integrals  Vsmr  sm  a,  Vp^  l)m  a,  V,,,,/  pmx,  etc.  are  given. 

Boolean  variable  noSpin  is  true  if  there  is  no  spin,  and  false  otherwise.  H*®,  is  given. 

if  noSpin  then 

Let  matrices  \H,j}  and  \S,j}  be  of  size  VortlValoms  x  VorbVa,oms. 
else 

Let  matrices  [Hu]  and  [S/y]  be  of  size  2VorbVatoms  x  2NotbNMoms. 

end  if 

Let  i  =  0. 

for  o'  in  s,  px, py, ...  do 
for  in'  —  I  to  Vatoms  do 

Increment  /  by  1  and  set  j  =  0. 

for  o  in  s,px,py, ...  do 
for  111  =  1  to  VaIorns  do 
Increment  j  by  I. 

Estimate  Ssw„m. 
if  in'  =  m  then 

Let  H,ymidm  =  E0z,„S,y5 

else 

Calculate  Vg^sm  as  a  linear  combination  of  Vslrf  sm  a,Vpmi  ,ml.o,  etc. 

Let  Hglmig„ |  =  2  ( EtTZ-m  EgZ'm)^o'm’dm  T  ^o'm'dnr 

end  if 

if  noSpin  then 

Let  I  =  i,  J  =  j ,  H/j  ~  H <7 ni  otn.  and  S/j  =  S(-/n/(-)nl . 

else 

for  S'  in  f,  J,  do 

for  5  in  1  do 

Let  /  =  2i  —  1  if  S'  =T,  and  let  /  =  2i  otherwise. 

Let  J  =  2j  —  I  if  S  =t,  and  let  J  =  2 j  otherwise. 

Let  o'  =  o' S'  and  o  =  dS. 

Let  Hu  =  H(ymi5mSSs'  +  H^n,om  and  StJ  =  S,y m8Ssi 

end  for 
end  for 
end  if 

end  for 
end  for 
end  for 
end  for 

Solve  generalized  eigenvalue  problem  \HU] C'  =  £,[S/y]C'. 


Figure  16.  Tight-binding  method  for  a  finite  system. 
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Like  H0im>om,  H0'm'0mk  can  be  broken  down  into  tenns  containing  two-center  integrals  such  as 
vsmr,sm,a,  vPm' ,Pm,a ,  vPm' ,Pm,m  and  VSm' ,Pm,a  ■  One  begins  by  expressing  it  explicitly  as  a  double 
Bloch  sum  ( 54): 


H 


o'm'omk 


=  lim 


wcell  wcell 

1  V-  V-  aik.(RJ.  +  dm-R.,  — dm,) 


Wcen-^00  Aicen  L  ■ 

j= i  j  =  i 


and 


X  /  "  R;'  “  dm')Hie-^ozm{r-Rj  -  dm)d3r , 

GO 

=  ^  eik  (Rt+dm~V)  |  <t>0'ZmX r  -  dm')Hle-(poZm(r  -  R;  -  dm)d3r . 


(89) 


The  integral  in  equation  89  is  essentially  identical  to  the  integral  that  defines  H0imiom  in 
equation  55,  with  Rm/  =  dm>  and  Rm  =  R;  +  dm,  and  the  process  for  decomposing  this  integral 
into  two-  and  three-center  integrals  is  the  same  as  that  for  H0rmr0m.  As  with  finite  systems,  all 
of  the  three-center  integrals  are  neglected.  Most  of  the  two-center  integrals  are  as  well, 
especially  the  ones  for  which  |R;  +  dm  —  d m>  |  »  0.  This  reduces  the  infinite  summation  in 
H o'm'omk  to  a  finite  one.  For  example,  in  a  crystal  with  the  zincblende  structure,  the  vectors 
connecting  atom  m  —  1  to  its  four  nearest  neighbors  are 

Ri  +  dx  —  d2  =  a(+ex  +  ey  +  ez)/4  , 

R2  +  dx  -  d2  =  a(+ex  -  ey  -  ez)/4, 

R3  +  dx  -  d2  =  a(-ex  +  ey  -  ez)/4, 
and 

R4  +  d1  -  d2  =  a(-ex  -  ey  +  ez)/4,  (90) 

where  a  is  the  lattice  constant  of  the  crystal.  If  Ssls2  =  0,  then  if  k  =  kxex  +  kyey  +  kzez  and 
all  but  the  nearest  neighbors  are  neglected  (55,  72): 

}-[  ,  —  TZ  |"gt(fca:+  ky+kz)a/4  p  gi(kx— /cy— kz)a/4  i  gi(— kx+  ky— kz)a/4 

SiSZK  51,5Z,(7L  /a  i  \ 

_(_  gi(— kx—  ky+kz)a/ 4] 


Also,  if  only  nearest  neighbors  are  taken  into  account,  diagonal  elements  of  H0rmr0mk  are  the 
atomic  eigenvalues  for  the  orbital  on  atom  m,  i.e.,  Homomk  =  E ozm  ■  The  overlap  matrix 

s o'm'omk  is  detennined  much  as  H0'm'omk »  he., 
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1 


and 


^cell  ^cell 


S  ,  ,  i,  =  lim  - V  V  eik  (R'+dm  V  dm') 

;=i  j'=i 


x  0o'z  ,(r  -  -  dm')0ozm(r  -  R;  -  dm)d3r 

J  771 


00 

=  £  eik  <5/+d”-dn.’)  J  tfyZm,(r  -  dmO0ozm(r  -  R;  -  dm)d3r . 

7  =  1 


(92) 


The  integral  in  equation  92  is  essentially  identical  to  the  integral  that  defines  S0imi0m  in  equation 
56,  with  R mt  =  d m>  and  Rm  =  R;  +  dm.  Figure  17  shows  a  pseudocode  implementation  of  the 
tight-binding  method  for  an  infinite  crystal. 
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The  spin-orbitals  s  "f,  s  [  px  f,  px  etc.  for  the  basis  set  are  given.  The  number  of 
spin-orbitals  in  the  set  is  Norb. 

The  number  of  atoms  per  unit  cell  Axioms  is  given. 

Values  of  two-center  integrals  Vjm- „„.CT,  V,, .CT,  Vpm>  pm  „,  etc.  are  given. 
Wavevectors  k| ,  ki, . . . ,  k,vk  are  given. 

for  k  in  k|  .k2,...,k/vk  do 

Let  matrices  [Hu]  and  [S/j]  be  of  size  NorbNMQm%  x  NOTbNMOm%- 
Let  /  =  0. 

for  o'  in  s  f,  s  [,  px  px  j, . . .  do 
for  m'  =  1  to  A/atoms  do 

Increment  /  by  1  and  set  7  =  0. 

for  o  in  s  ],  s  px  f,  px  J., . . .  do 

for  m  =  1  to  Axioms  do 
Increment  J  by  1 . 

Let  H,jmtomV  =  0  and  =  0. 

for  R j  in  R1(  Rj, ...,  Rmax  do  {where  |R||  <  |R2 1  <  <  Rmax} 

Let  atom  M  be  the  atom  in  the  crystal  located  at  R j  -(-  d,„. 
Determine  H(,imi,M  and  S0i„,i0m  as  one  would  for  a  finite  system. 

Add  eik  (fty+d--d™')f/oWoM  to  H^om k. 

Add  eik  < Rr +dm ~d<»' ) S0<m'oM  to  SoWom k. 

end  for 

Let  Hu  =  HoWomk,  and  Su  = 

end  for 
end  for 
end  for 
end  for 

Solve  generalized  eigenvalue  problem  [Hu]C'k  —  £7k[S/./]C'k. 

end  for 


Figure  17.  Tight-binding  method  for  an  infinite  crystal. 


Here  are  examples  of  how  the  special  two-center  integrals  Vsm>  smcT,  Vpm'.pm.o’  etc.,  maY  be 
estimated  through  closed-form  empirical  fonnulas.  One  such  formula  ( 56)  is 


V, 


'  R 


(0)  N’VoA 


o  m  ,om,X 


K(0) 

o'm' ,om,X  ’ 


(93) 


where  o  E  {s,  p, ... },  A  E  {a,  n),  R^m'  kS  the  distance  between  atoms  m'  and  m  in  the  absence  of 
strain,  gmA  is  the  strain- free  value  of  that  is,  the  value  of  V0'm'  om^  evaluated  at 
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anC*  Vo  i0A  is  an  empirical  parameter  that  is  typically  in  the  range  of  1  to  4.  l/^,  om  X  is  an 
empirical  parameter  itself.  Another  formula,  from  Mehl  and  Papaconstantopoulos  (19),  is 


V, 


o'  m!  ,om,A 


2 

=  (e  i  i  +  f  i  iR  i  +  f  i  iR2  (R  A 

Vco  oA  ‘  J  o  oA/v mm  '  J  o  oAlYmm 1  c\lxmm  J 


(94) 


where  e0toX,  f  0'0x,  f  0'0A’  anc'  g0’oA  are  fitting  parameters,  and  Fc(Rmrn /)  is  a  cutoff  function, 

1 


Fc(Rmm 0  = 


1  +  exp  [(flmm'  Rmm’)/Ij\ 

Mehl  and  Papaconstantopoulos  typically  take  L  to  be  half  the  Bohr  radius. 


(95) 


Off-diagonal  elements  of  the  overlap  matrix  may  be  determined  through  closed-fonn  empirical 
formulas  as  well.  For  example,  Papaconstantopoulos  et  al.  ( 67)  formulate  this  dependence  as 
follows: 


S 


o'm'  ,om,A 


o' o  Q-o' oARmm'  "h  b0'0AR  mm'  "f"  ^ o' oA^mm'^) 

x  e-d2o'oARmm'Fc(Rmm,)  , 


(96) 


where  ci(/oA.  b0ioA,  c0<oX,  and  d(/oA  are  fitting  parameters.  The  Kronecker  delta  80/0  ensures 
that  the  preceding  expression  is  consistent  with  the  condition  S0imom  =  80r0  (where  m'  =  m), 
which  the  overlap  matrix  must  satisfy. 


Such  formulas  are  a  part  of  how  the  tight-binding  method  takes  strain  into  account,  since  they 
approximately  account  for  changes  in  interatomic  distance  due  to  strain.  Strain  also  leads  to 
changes  in  the  orientations  of  atoms  relative  to  one  another  (i.e.,  the  directions  of  Rmm/  for  pairs 
of  atoms  m'  and  m),  and  this  is  taken  into  account  through  the  direction  cosines  lR,  mR,  and  nR 
used  in  the  decomposition  of  V0'm'0m  into  linear  combinations  of  Vsm'  sma,  VSm',pm,o >  etc-  (and 
also  into  any  decomposition  of  S0imi0m  into  linear  combinations  of  Ssrn>  pm  a,  Spm>  prn  (T,  etc.). 

It  can  be  shown  that  the  diagonal  elements  of  H0imr0rn  or  H0im'orrik  also  depend  on  the  strain. 
According  to  Boykin  et  al.  (33),  if  4>0zm  (r  —  Rm)  is  a  Lowdin  orbital  as  defined  in  equation  57, 
then 


J  OZy, 


—  FoZm  So'  Tim'  C0 


Ew  (0) 

coZm+c0'z  , 
m 


(97) 


where  the  superscript  “(0)”  again  denotes  the  strain-free  version  of  a  quantity,  and  Co'm'om 's  an 
empirical  parameter.  This  is  not  the  only  scheme  for  accounting  for  the  effects  of  strain  on 
diagonal  elements.  For  example,  Mehl  and  Papaconstantopoulos  (19)  take  the  diagonal  elements 
to  be 
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Homom  & oZm  ”1”  boZmPm  ~f  ^ oZmPm  ~f  ^ oZmPm  ■> 


(98) 


where  aoZm,  boZjn,  coZjn,  and  doZm  are  fitting  parameters,  and 


Pm  = 


^  '  exP  (  1 

m' 


(99) 


The  summation  with  index  m!  is  over  the  neighbors  of  atom  m,  and  ^zmz  ,  is  a  fitting 

parameter.  (Mehl  and  Papaconstantopoulos  (19)  fit  their  parameters  to  results  from  ab  initio 
calculations.)  Other  methods  for  accounting  for  strain  on  diagonal  elements  have  been  discussed 
by  Jancu  and  Voisin  (73)  and  Niquet  et  al.  (74). 


The  Slater-Koster  tight-binding  method  has  been  used  on  systems  with  about  a  million  atoms. 

For  example,  it  has  been  used  on  a  dome-shaped  Ino.6Gao.4As  quantum  dot  embedded  in  a  GaAs 
matrix.  The  dot  was  30  nm  in  diameter  and  5  nm  in  height,  and  the  simulation  domain  was 
40  x  40  x  15  nm.  The  dot  itself  contained  about  718,000  atoms,  with  the  GaAs  matrix 
constituting  the  remainder  of  the  million  atoms  in  the  simulation.  The  distribution  of  the  In  and 
Ga  cations  within  the  dot  was  random.  For  each  cation  (Ga  or  In)  and  anion  (As),  the  set  of  basis 
functions  contained  three  p-orbitals  (px,  py,  and  pz)  and  two  orbitals  with  s-symmetry,  denoted 

s  and  s* .  Since  spin  was  taken  into  account,  the  number  of  basis  functions  per  atom  doubled 
from  5  to  10.  For  a  given  distribution  of  Ga  and  In  atoms  within  the  dot,  detennination  of  the 
energy  eigenvalues  of  the  system  took  about  25  min  on  3 1  processors  of  a  cluster  of  Pentium  III 
933-MHz  CPUs.  If  the  corresponding  eigenvectors  had  been  computed,  the  computation  time 
would  have  doubled  (56).  As  seen  from  the  pseudocode  of  the  tight-binding  implementations 
shown  in  figures  5  and  6,  the  computational  resources  needed  for  the  method  scales  with  the 
product  of  the  number  of  orbitals  in  the  basis  set  and  the  number  of  atoms  in  the  system. 

Whether  it  scales  linearly  or  as  0(N 3)  depends  on  the  choice  of  eigenvalue  solver  (21,  46). 


3.3  How  Non-Self-Consistency  Affects  How  Strain  Is  Taken  Into  Account 


The  atomistic  methods  previously  described  are  not  self-consistent,  that  is,  their  approximations 
of  the  effective  potential  of  the  one-electron  Schrodinger  equation  do  not  depend  on  the  charge 
density  or  the  one-electron  wavefunctions.  Computationally,  this  is  a  great  advantage,  since  it 
means  that  the  Schrodinger  equation  does  not  have  to  be  solved  iteratively.  However,  as  will  be 
explained  shortly,  this  entails  a  restriction  on  the  physics  that  the  effective  potential  Text, eft  maY 
take  into  account.  For  a  system  of  M  atoms  (where  M  may  be  infinite),  Fext,eff(r-  (R)m)  may  be 
written  generally  as 


^ext,eff(U  (R) M ) 


M 

1  y  Ztq2 
47re0Zj|r  — RJ 

i=i 


+ 


P(r') 

|r-r'| 


fi  I"  T  Pother  j 


(100) 


where  the  operator  Fothcr  is  determined  by  one’s  choice  of  ab  initio  method.  This  can  be 
contrasted  with  the  decompositions  of  Fext,eff  into  a  sum  of  per-atom  contributions  in  equation  3 1 
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for  the  empirical  pseudopotential  method  and  equations  62  and  63  for  the  Slater-Koster 
tight-binding  method.  Despite  the  contrasting  details  in  the  earlier  descriptions  of  these  two 
methods,  we  remark  that  they  both  follow  similar  algorithmic  and  fundamental  theoretical 
structures.  Namely,  in  the  generalized  effective  potential  in  equation  100,  only  the  first  term,  the 
contribution  from  the  M  nuclei  of  the  system,  naturally  decomposes  into  a  sum  of  per-atom 
contributions.  One  could  decompose  the  second  and  third  tenns  in  Vextcff  by  assuming  the 
following: 


and 


such  that 


M 


P(r)  =  ^Pi(r)  , 


i=i 


V, 


other 


Zf=1  Vi 


t.other  ■ 


(101) 
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bext,eff(D  Wm) 
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[.other 


(103) 


Barring  an  accounting  of  some  nonlocal  effects  such  as  spin-orbit  coupling  (e.g.,  Chelikowsky 
and  Cohen,  [35]),  this  is  largely  what  the  atomistic  methods  described  are  effectively  doing. 
Furthermore,  in  these  methods,  p,  (r)  is  generally  dependent  on  the  species  of  atom  i  (i.e., 
whether  it  is  Si,  Ga,  nitrogen  [N],  etc.)  rather  than  its  location  in  the  system  Rj.  This  is  a 
problem,  for  example,  in  piezoelectric  systems,  where  an  electric  field  due  to  strain  may  lead  to 
the  transfer  of  charge  from  atoms  near  one  region  of  the  system  to  atoms  near  another  region, 
even  if  the  atoms  in  each  region  are  of  the  same  types. 


This  problem  can  be  overcome  by  treating  this  field  as  if  it  were  external  and  adding  the 
potential  due  to  this  field,  qVstat(r,  e),  where  e  is  the  strain  if  the  field  is  induced  from 
piezoelectricity,  to  the  effective  potential  (34,  75-78),  that  is, 


M 


^ext,eff(T;  {R}m)  —  /  ^  ^U6 


Ztq2 


f  Pi  (r') 
J  |r  —  r' 


.  ,  „  Jr  -  Ril 

i=i  L 

M 

+  ^  Mother  +  qVstat(r,  e). 
[=1 


dV 


(104) 


Fstat(r,  e)  can  be  determined  from  a  continuum  mechanical  calculation  (75,  79,  80),  provided 
that  the  strain  e  in  the  system  has  already  been  determined  from  the  changes  in  atomic  positions. 
Possible  methods  to  detennine  this  strain  are  in  appendix  B  of  this  report. 
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4.  The  k  •  p  and  Envelope  Function  Methods 


While  the  atomistic  empirical  methods  are  far  less  computationally  expensive  than  the  ab  initio 
methods,  their  calculation  time  still  scales  with  the  number  of  atoms  in  the  simulated  system. 
This  is  not  an  issue  for  the  k  ■  p  and  envelope  function  methods,  which  do  not  take  atomic 
positions  as  input  at  all.  The  k  ■  p  method  (81)  was  originally  developed  to  estimate  the 
electronic  band  structure  of  bulk  periodic  semiconductor  crystals,  especially  the  band  structure 
near  the  T  point.  The  envelope  function  approximation,  however,  extends  the  method  to  systems 
that  are  not  periodic  in  all  dimensions,  such  as  quantum  wells,  wires,  and  dots.  In  this  extension, 
the  effective  potential  within  a  given  material  is  still  assumed  to  behave  approximately  as  the 
effective  potential  would  in  the  bulk  form  of  the  material.  Instead  of  taking  strain  into  account 
through  atomic  positions,  it  is  taken  into  account  through  an  operator  added  to  the  effective 
Hamiltonian,  an  operator  that  takes  the  small  strain  tensor  from  continuum  mechanics  as  input. 
These  approximations  mean  that  the  envelope  function  approximation  is  essentially  a  continuum 
approach  rather  than  an  atomistic  one. 

4,1  General  Formulation  of  the  k  ■  p  Method  for  Bulk  Crystal 

The  k  ■  p  method  for  bulk  crystal  employs  the  Bloch  theorem,  so,  accordingly,  the  one-electron 
wavefunction  is  expressed  as  i/'tkC1')  =  elk'riqk(r),  where  u(k(r)  has  the  same  periodicity  as  the 
crystal.  Substituting  this  expression  for  xpik( r)  into  the  one-electron  Schrodinger,  equation  8 
yields  (81) 


+  Pext,eff(r-  Woo)  +  ^-k  •  p  +  uik(  r)  =  Eikuik(  r) 


(105) 


where  k  2  =  |k|2  ,  and  (R}oo  is  the  set  of  coordinates  and  atomic  numbers  of  the  nuclei  of  the 
infinite  bulk  crystal.  Alternately,  equation  105  may  be  written  as 


Hklep-nik(r)  =  [tfle-( r, Woo)  +  ^-k  •  p  +  uik(r)  =  Eikuik( r) 


(106) 


where  Hle-  (r,  (R}oo)  is  the  one-electron  Hamiltonian  for  an  infinite  bulk  crystal.  For  small  k, 
the  operator  is  “close”  to  Hle-,  so  the  eigenstates  of  H^e-  may  be  expanded  in  tenns  of  the 
eigenstates  of  Hle-  at  k  =  0: 

W;k(r)  =  c^k)um0(r)  .  (107) 


These  eigenstates  are  orthogonal,  that  is, 
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fn-Um'o  (r)umo  (r)  d3  r  =  8r 


(108) 


where  flc  is  the  volume  of  the  primitive  unit  cell  of  the  crystal.  Alternatively,  iqk(r)  may  be 
expanded  not  directly  in  terms  of  the  eigenstates  at  the  o  point  but  rather  in  terms  of  normalized 
orthogonal  linear  combinations  of  these  eigenstates,  i.e., 

Mik(r)  =  EmCmk)tfm(r),  t/f( r)  =  ZmCinUmoir)  • 


(109) 


If  one  substitutes  equation  109  into  equation  106,  multiplies  the  latter  equation  by  (7^/(r),  and 
integrates  over  a  unit  cell  of  the  crystal,  then  one  may  obtain  a  matrix  equation 

V  lA'P  rfDl  ^/Aik)  _  p  p(tk) 

ZjTTL  nm'm  vtWj co J^m  ^  m' 
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tfm'm(Woo)  =  /n  t/^(r)//le-(r,{R}00)t/rn(r)d3r, 


Pm'm  =  fn  l/^(r)pt/rn(r)d3r. 


(110) 

(HI) 

(112) 

(113) 


Alternately,  h2k2/2me  may  be  moved  to  the  right-hand  side  of  the  matrix  equation,  so  that 

£  ffX«Rudik)  =  (%  -  •  (114) 


where 


=  R-P« 


h 


(115) 


t/j(r)  is  often  chosen  such  that  Hlc-  UL(r)  —  Ei0Ui(r),  so  that  H*®  —  Em08m'm  (32,  81-83). 


If  the  effect  of  spin-orbit  is  taken  into  account,  one  may  then  replace  p  in  equations  106  and 
1 13  with  ft,  where  (32,  82) 


=  P  +  3^  [®  x  ■ 

When  k  is  small,  kft«kp.  (83). 


(116) 
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4.2  Example  k  ■  p  Formulations  for  Bulk  Crystals 

In  principle,  the  summation  in  the  expansion  of  iqk(r)  in  equations  107  and  109  is  infinite  and, 
accordingly,  so  is  the  matrix  H  ,  (or  Hr).  In  practice,  Hm’m  is  transformed  into  a  finite 
k  ■  p  Hamiltonian  matrix,  usually  one  with  fairly  small  dimensions,  e.g.,  6  x  6  or  8  x  8.  Two 
examples  of  how  this  may  be  done,  the  Kane  and  Luttinger-Kohn  formulations,  are  shown  in  the 
following  equations.  In  the  formulation  by  Kane  (53),  the  transformation  is  done  simply  by 
truncating  the  expansion  of  iqk(r),  including  only  the  terms  corresponding  to  the  one-electron 
wavefunctions  for  the  three  valence  bands  and  the  lowest  conduction  band  of  a  semiconductor 
with  a  diamond  or  zincblende  crystal  structure.  When  spin  is  taken  into  account,  consideration 
of  these  four  bands  leads  to  an  eight-term  expansion  of  iqk( r) — four  terms  for  spin-up  and  four 
for  spin-down — and  thus  an  8  x  8  matrix.  In  the  Luttinger-Kohn  formulation  (32),  the  expansion 
of  iqk(r)  is  infinite  in  principle,  but  the  six  terms  pertaining  to  the  three  valence  bands  are 
assumed  to  be  the  dominant  ones  in  the  expansion,  and  the  infinite  matrix  is  accordingly 
transfonned  into  a  6  x  6  matrix  where  the  effects  of  the  nondominant  terms  in  the  expansion  of 
iqk(r)  are  indirectly  taken  into  account  through  material  constants  called  the  Luttinger 
parameters. 

Kane  (S3)  studied  the  band  structure  of  indium  antimonide  (InSb),  whose  electronic  band 
structure  resembles  the  schematic  shown  in  figure  18.  A  is  the  difference  in  energy  between  the 
maximum  of  the  split-off  valence  band  and  the  common  maximum  of  the  other  two  valence 
bands,  labeled  as  “heavy  hole”  and  “light  hole.”  Those  labels  refer  to  the  effective  masses 
(discussed  in  detail  in  section  4.4)  of  the  electrons  or  holes  with  energies  in  the  corresponding 
bands.  Holes  are  vacancies  left  behind  when  an  electron  is  promoted  from  a  valence  to  the 
conduction  band,  and  they  act  like  positive  charge  carriers.  The  expansion  in  equation  109  is 
truncated  to  include  only  the  eight  states  at  k  =  0  corresponding  to  the  bands  shown  in  the 
figure.  The  two  states  corresponding  to  the  conduction  band  minimum,  one  for  spin-down  and 
one  for  spin-up,  have  the  same  radial  symmetry  as  atomic  s-orbitals,  and  these  states  are 

14 (r)  =  iS(r)  i  ,  and  l/5(r)  =  iS(r)  T  ,  (11V) 

where 

t=[J],andl=[J],  (118) 

and  S  (r)  is  spherically  symmetric  and  thus  depends  only  on  the  magnitude  of  r.  Accordingly,  it 
has  the  following  parity  properties: 

5((r1,r2,r3))  =  S((|r1|,  |r2|,  |r3|))  .  (119) 
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Figure  18.  Schematic  band  structure  near  k  =  0  of  a  typical  semiconductor  with 
a  diamond  or  zincblende  crystal  structure.  Eg  is  the  band  gap  energy. 
A  is  the  difference  in  energy  between  the  maximum  of  the  split-off 
band  and  the  common  maximum  of  heavy  and  light  hole  bands.  The 
left  and  right  halves  of  the  horizontal  axis  indicate  the  magnitude  of 
k-values  pointing  along  certain  crystal  directions.  (In  this  qualitative 
schematic,  the  actual  directions  are  not  important.)  While  the  heavy 
hole,  light  hole,  and  split-off  bands  are  typical  for  such  a 
semiconductor  and  generally  have  their  maxima  at  k  =  0  as  shown, 
the  actual  minimum  of  the  conduction  band  may  be  different  from 
what  is  shown  in  this  schematic  (77). 


The  six  valence  states  at  k  =  0,  which  have  the  same  symmetries  as  atomic  p-orbitals,  are 


U2(  r) 


X(r)-iHr)  T 
V2 


u6( r)  =  - 


X(r)+iF(r) 

vf 


i. 


(120) 


U7(r)  =  Z(r)  T  ,  (121) 

t/8(  (122) 

where  the  functions  X,  Y,  and  Z  have  the  following  parity  properties: 

X((ri,r2,r3))  =  sgn(r1)X((|r1|,  \r2\,  |r3|))  ,  (123) 

T((ri,r2,r3))  =  sgn(r2)T((|r1Ur2Ur3|))  ,  (124) 

and 

Z((ri,r2,r3))  =  sgn(r3)Z((|r1|,  |r2|,  |r3|))  ,  (125) 

where  sgn(x)  is  the  sign  of  x,  i.e.,  sgn(+x)  =  +1.  Kane  assumed  that  k  •  ft  w  k  •  p,  since  k  is 
taken  to  be  small.  Kane  initially  restricted  the  Bloch  wave  vector  k  to  points  along  the  z- 
direction  (i.e.,  the  r3  axis)  so  that  the  matrix  would  simplify  to  the  following  block 
diagonal  format: 


t/3(r)  =  Z(r)  i, 


ir  _  *(r)+iK(r)  ^ 
~  V2  T’ 
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where 
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and 


Ep=l  Xl(r)4-(r,{R}00)I(r)d3r, 

J  nc 

=  f  Zt(r)Hle- (r,  {R}DO)T(r)d3r  , 
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PKane  =  f  Sf(r)  p3Z( r)d3r, 

me  Jnc 

(130) 

3hi  f  j.  ,  .  9 kextgff  d^ext.eff  „  ,3 

4m|^Jfl/t(r)[  an  an  p'P(r)dr- 

(131) 

where  p[  is  a  component  of  vector  operator  p.  PKane  is  called  Kane’s  parameter  (32).  Here, 
HKane,4x4  is  written  under  the  assumption  that  equation  114,  the  k  ■  p  matrix  equation  with  the 
term  —  h2k2 /2me  on  the  right  side,  is  being  used.  If  k  points  in  a  general  direction,  a 
transformation  can  be  applied  to  rotate  the  original  coordinate  system  so  that  k  points  along  the 
z-direction,  which  makes  HKane  block  diagonal.  Since  S(r)  is  spherically  symmetric,  it  is 
unaffected  by  this  rotation.  The  rotation  causes  Z(r),  T(r),  Z(r),  T,  and  -l  to  be  replaced  by 
Z'(r),  T'(r),  Z'( r),  T',  and  -1'  in  equations  1 17,  120-122,  128,  and  131.  If  k  is  expressed  in  the 
original  coordinate  system  in  spherical  coordinates  ( k ,  6,  0),  the  primed  functions  are,  in  tenns 
of  the  original  Z(r),  Y  (r) ,  and  Z(r), 


Z'(r)  =  X(r)  cos  6  cos  <p  +  Y( r)  cos  0  sin  0  —  Z(r)  sin  0  , 


Y'( r)  =  —  Z(r)  sin  <p  +  Y( r)  cos  0 , 
and 

Z'(r)  =  Z(r)  sin  6  cos  0  +  T(r)  sin  6  sin  0  +  Z(r)  cos  6  .  (132) 


The  primed  spinors  are 
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and 


e  1(^/2cos(0/2) 
evp/2  s[n(g/2) 


— e  ^Z2  sin(0/2) 

e\4>/2  cos(q/2) 


H Kane, 4x4 >s  itself  in  a  block  diagonal  format,  with  one  block  containing  the  first  three  rows  and 
columns  and  the  other  block  being  a  1  x  1  “matrix”  with  Ep  +  A°/3  as  its  sole  element.  Thus, 
one  may  define  two  eigenproblems,  one  for  each  block,  as  follows: 
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(134) 


(135) 


(£„  +  4”/3 )C«k)  =  (Eik  -  |£)  C«k)  .  (136) 

In  order  for  the  highest  valence  band  extremum  to  be  zero,  Ep  is  set  to  —  A°/3.  This  leads  to  the 
following  solution  for  Eik: 


Eik  ~  Ef iji  k  — 


h 2  , 

9 - k 

2  me 


(137) 


This  is  the  heavy  hole  band  according  to  the  Kane  formulation.  Note  that  this  band  curves 
upward  as  a  function  of  k,  which  contradicts  experiment  (32,  83)  (and  the  schematic  band 
structure  in  figure  18).  This  is  a  consequence  of  truncating  the  expansion  of  uik(r).  Because  of 
spin,  this  energy  band  is  doubly  degenerate,  and  the  periodic  parts  of  the  two  one-electron 
wavefunction  that  corresponds  to  this  band  are,  for  k  pointing  along  a  general  direction, 
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uhhik(x)  =  t/4(r) 


X'(r)+iy'(r)  v 
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(138) 


where  the  primed  quantities  are  defined  in  equations  132  and  133. 


The  eigenvalues  of  HKane  3x3  may  be  determined  numerically.  However,  some  insight  into  the 
physical  meanings  of  Es  and  A0  may  be  obtained  by  estimating  closed-form  expressions  from 
these  eigenvalues,  similar  to  the  one  for  Ehh.k-  The  characteristic  polynomial  of  HKane  3x3  is 

ihkOhk  -  EsXEq k  +  A°)  -  k2P2(E(k  -  2A°/3)  =  0  ,  (139) 
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where  Eik  —  Eik  —  h2k2 /2me,  and  the  equation  Ep  —  —  A°/3  has  already  been  substituted  into 
H Kane, 3x3-  If  k  =  0,  the  solutions  of  the  characteristic  polynomial  are  E'i0  =  Es,  0,  — A°.  These 
solutions  are  presumed  to  be  extrema  of  the  energy  bands.  For  simplicity,  these  bands  are  also 
presumed  to  be  parabolic  in  the  neighborhood  of  k  —  0.  Accordingly,  trial  solutions  for  this 
polynomial  for  k  A  0  are  taken  to  be  of  the  form  E'ik  —  Es  +  e^k2, 0  +  e^k2,  — A°  +  e^/c2,  where 
et  is  a  parameter  to  be  determined  (32).  If  these  trial  solutions  are  substituted  into  the 
characteristic  polynomial  and,  then,  any  resulting  terms  containing  powers  of  k  higher  than  2  are 
ignored,  then  for  i  —  c,  lh,so, 

_  Planers  +  2A°/3)  2 P2ane  ^Kane  / 1  a 

6c  "  ES(ES  +  A°)  ’  6lh  3ES  ’  gso  ~  3  (Es  +  A°)  ’  1  ' 

and  the  energy  bands  are  approximately  (32,  83) 
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where  Eck,  Em  k,  and  Eso  k  are  the  energies  for  a  given  k  in  the  conduction,  light  hole,  and  split- 
off  bands,  respectively.  By  comparing  these  results  to  the  qualitative  band  structure  in  figure  18, 
one  can  confirm  that  Es  —  Eg  and  A0  =  A.  One  may  also  note  from  equation  131  that  if  the 
spin-orbit  coupling  were  nonexistent,  parameter  A  would  not  exist  either.  Eg,  A,  and  PKane  maY 
be  detennined  empirically,  and  a  compilation  of  these  parameters,  along  with  some  discussion  as 
to  how  they  have  been  obtained,  may  be  found  in  a  survey  by  Vurgaftman  et  al.  (84). 


Because  of  spin,  the  energy  bands  Eck,  Eih  k,  and  Eso  k  are  doubly  degenerate.  For  k  pointing 
along  a  general  direction,  the  periodic  parts  of  each  one-electron  wavefunction  iqk(r)  are 
determined  to  be 


%k(r)  =  ajU1  (r)  +  bjU2( r)  +  CjU3( r) , 


=  ajiS{r)  l'+b,x'W^r'M1+c,Z'(r)  V  , 


Uj2k(r)  =  ajU5( r)  +  bjU6( r)  +  CjU7( r)  , 
and 

a;iS(r)  V-  bj  i-+  CjZ'(r)  T' , 


(142) 
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where  (a,-,  bj,  Cj)  is  an  eigenvector  of  HKane  3x3,  and  j  —  c,  Ih,  so.  When  k  =  (0,0,0),  then 

ac  —  1,  bc  —  0,  cc  =  0  , 

am  =  0?  blh  —  Vl/3 ,  Cih  =  -J2/3  , 

^SO  —  hso  —  V2/3  ,  Cso  —  —yJ\/3 


Pseudocode  of  the  Kane  formulation  of  the  k  ■  p  method  is  shown  in  figure  19. 


Parameters  Eg,  A  and  fkanc  are  given. 

Wavevectors  k|,  k2, . . .,  k,\k  are  given,  where  k,  is  reasonably  close  to  T  point. 
H Kane. 3 x 3 {k,  Es ,  Ep ,  A,  /kanc )  is  given. 


for  k  in  k|  .k2,...,kwk  do 

Let  k  =  |k|.  {k  in  spherical  coordinates  is  (k,8,<j>).} 


LetX'(r)  =  X  (r)  cos  0  cos  0-FF(r)  cost?  sin  0  —  Z(r)sin0. 

Let  t"(r)  =  —  X(r)sin0  +L(r)cos0. 

Let  Z'(r)  =  X(r)sin0cos0  -t-Z(r)sin0sin0  -Z(r)cos0. 

(A(r)  has  spherical  symmetry  and  isn’t  changed  by  a  coordinate  rotation.} 


Let  f= 


Let  i'= 


e  '^,2cos(0/2) 
ei0/2sin(0/2) 

— e“'^/2sin(0/2) 
e1^/2  cos(0/2) 


{Defining  spin-up  in  spherical  coordinates.} 

.  {Defining  spin-down  in  spherical  coordinates.} 


{For  the  conduction,  light  hole  and  split-off  bands:} 

Let  H  =  HKane,3x3(*)£g!_A/3,A,flcane)- 
Find  eigenvalues  E'ik  and  eigenvectors  C,k  of  H. 

Let  a,  =  Cf,  b,  =  Ctk.  and  c,  =  Cf. 

for  i  in  c,lh,so  do 

Let  £,k  =  E'ik  +  trk1  j2mc. 

Let  w;ik(r)  =  a,iS(r)  {'  +fc,{[X'(r)  -  iZ'(r)]/v/2}  T'  +c,Z'(r)  {'. 
Let  «,2k(r)  =  «,iS(r)  {'  -fc/{[X'(r)  +  \Y'(r)\/V2}  }'  +c,Z'(r)  }'• 

end  for 

{For  the  heavy  hole  band:} 

EhhV.  =  h2kr  /2me. 

«/,/,ik(r)  =  {[X'(r)  +  iZ'(r)]/^}  }'• 

«AMk(r)  =  {[X'(r)-iZ/(r)]/v^}  {'. 

end  for 


Figure  19.  Pseudocode  for  Kane  formulation  of  the  k  ■  p  method. 


(143) 
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The  determination  of  the  Hamiltonian  matrix  HKane  does  not  require  a  detailed  knowledge  of  S, 
X,  Y,  or  Z;  rather,  only  the  parity  properties  shown  in  equations  1 19  and  123-125  are  required. 
This  causes  most  of  the  calculated  matrix  elements  to  be  zero.  The  integrals  that  are  nonzero 
become  functions  of  empirical  parameters.  This  is  a  common  feature  of  k  ■  p  formulations. 

While  Kane  (53)  simply  truncated  the  expansion  in  equation  107  to  obtain  a  finite  matrix  for  his 
k  ■  p  formulation,  Luttinger  and  Kohn  (82)  took  a  different  approach  that,  in  principle,  accounted 
for  all  the  terms  in  the  expansion.  They  broke  the  expansion  into  two  sums: 

A  B 

cmk)um o(r)  +  ^  cfk)ui0(r)  .  (144) 

l 

The  first  summation  is  over  the  class  A  eigenstates;  that  is,  the  few  states  whose  eigenenergies 
Em o  are  near  the  bottom  of  the  band  gap.  The  second  summation  is  over  the  class  B  eigenstates, 
which  are  the  remaining  states.  Here,  Um  =  um0,  so  equations  107  and  109  are  equivalent  and 
c£k)  =  C£k) .  The  class  A  states  may  be  expressed  as  (32) 


W;k(r)  =  ^ 


and 


(145) 


X,  Y,  and  Z  have  the  same  meanings  they  do  in  the  fonnulation  by  Kane  (53)  previously 
discussed.  These  states  are  also  the  periodic  parts  of  the  one-electron  wavefunctions 
corresponding  to  the  heavy  hole,  light  hole,  and  split-off  bands  at  k  =  0,  just  as  in  the  Kane 
model.  (See  equations  138,  142,  and  143  for  comparison.)  They  have  the  same  symmetry  as  the 
p-eigenstates  of  an  atom  where  spin-orbit  coupling  is  taken  into  account  (35,  pp.  99,  197).  The 
states  in  equation  145  all  belong  to  the  valence  bands.  In  the  Luttinger-Kohn  formulation,  all 
conduction  band  states  are  within  class  B .  Because  of  the  parities  of  X,  Y,  and  Z,  pm'm  =  0  for 
the  class  A  states. 
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Equation  1 14  is  rewritten  in  matrix  form  as 


(H°  +  H1)C(ik)  =  Et 


ik 


h2k2' 
2  meJ 


C(ik) 


where  the  elements  of  H°  and  H1  are  such  that 


H m'm  E-mO^rn'm  > 


and 


T, 1 

nm'm 


k  P  r 


Luttinger  and  Kohn  (82),  here,  apply  a  canonical  transformation,  where 

c(tk)  _  eSB(ik)5 

and  then 


e“s(H°  +  H1)esB(ik)  =  Ei} 


h2k2' 
2  mP 


ik 


B®, 


A  new  matrix,  Hreduced,  is  defined  as 


H 


reduced 


h2k2 


e~s(H°  +  H1)es  + - 1 , 

2  mP 


where  I  is  the  identity  matrix  and 

jj reduced g(tk)  _  £-.kg(ik)_ 

The  first  tenn  of  Hreduced  may  be  expanded  as 

e~s(H°  +  H1)es  =  H°  +  H1  +  [H°,  S]  +  [H1,  S] 

+  ^  [[H°,  S],  S]  +  ^  [[H°,  S],  S]  H — , 

where  the  commutator  [A1;  A2]  =  A  ;  A2  —  A2AX.  The  matrix  S  is  chosen  so  that 

H1  +  [H°,  S]  =  0  , 
and 


(146) 

(147) 

(148) 

(149) 

(150) 

(151) 

(152) 

(153) 


Hreduced  ^  H0+i[Hl  S]  +g!,  (154) 

This  removes  coupling  terms  between  classes  A  and  B  that  are  first  order  in  k.  Given  the 
definitions  of  H°  and  H1,  equation  153  becomes 


k  •  Pm'm  3"  (.Em1 0  ^mo)^m'm  ~  0  >  (155) 

me 

and  the  elements  0f  Hreduced  become  (81),  for  indices  m' ,  m  ranging  from  1  to  6,  the  index 
values  pertaining  to  class  A, 
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irfeduced 

nm'm 


(k-Pm';)(k-P;?n)  |  (k-pm>t)(k-ptm)' 
^mo  —  Eio  Em>  o  —  El0 


(156) 


The  summation  index  l  only  ranges  over  the  class  B  states  because  pm'm  vanishes  for  states  in 
class  A.  If  it  is  assumed  that  the  states  in  class  A  all  share  at  least  approximately  the  same 
eigenvalue,  i.e.,  Em0  ~  Em'0  =  E0,  then  the  previous  equation  simplifies  to  (32) 


jjLK  (r  ,h2k2\r  ,  ^  V  ^k'Pmh)(k-pim) 

m'm  ~  v  m°  2m  J  m’m  2me  L  E„-Em 


(157) 


where  H~,/  is  called  a  Luttinger-Kohn  Hamiltonian  matrix,  and 


U£k(r)  ~  ^5mk)umo(0  ■ 


(158) 


The  summation  in  equation  157  containing  matrix  elements  involving  class  B  states  is  expressed 
in  terms  of  the  empirically  detennined  Luttinger  parameters  Yi,  ¥2,  and  Y?,-  These  can  be 
expressed  as 

Ki  =  “f^t4"+2Bo)> 


where 


<A0-Bo) 


_ _  r 

3h2  0  ’ 


(159) 


A  ft2  ft2  V  P^Pn_ 

°  2  me  m2  Z_i  E0  —  El0 

R  _  h2  |  ftZ  yn  PiM? 

0  2 me  mlLl  E0-El0  ’ 


'-I 


_(1)_(2)  ,  (2)  Cl) 

X  Pll  Pl2  P21  Pll 


(160) 


Eq  -  El 


and  Pill  is  the  ith  component  of  the  vector  p lm. 


Chuang  (32)  has  expressed  a  Luttinger-Kohn  Hamiltonian  matrix  as  follows: 
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(161) 


H 


LK 


where 


and 


P  +  Q 

-5 

R 

0 

V2P 

-5+ 

P-Q 

0 

R 

-V2<2 

77725 

P  + 

0 

P-Q 

S 

73725+ 

V2Q 

0 

p+ 

s+ 

P  +  Q 

-V2P  + 

—  5+/V2 

-5+/V 2 

-V2Q+ 

JV2S 

-yf2R 

P-A 

0 

V2P+ 

^/y2s^ 

-5/V2 

0 

P-A 

p  =  - 

£<«** 

:2  +  /Cl)  , 

Q  =  - 

h2y2  „  , 

(/q  +  k 

2m  e 

I  -  2k2)  . 

R  =  - 

h2 

ikl-k2) 

1  +  i2v/3y3 

kik2)  , 

S  =  - 

ft2r3V3(fel- 

~  ik2)k3  . 

m„ 


(162) 


Unlike  the  Kane  formulation,  the  Luttinger-Kohn  k  ■  p  formulation  is  written  for  Bloch  wave 
vectors  k  of  arbitrary  direction,  so  the  pseudocode  for  it,  shown  in  figure  20,  is  trivially  simple. 


Luttinger  parameters  ft ,  ft,  73  are  given. 

A  is  given. 

Wavevectors  ki ,  k?, . . . ,  k,vk  are  given. 

Luttinger-Kohn  Hamiltonian  matrix  H(k,  ft  ,ft,  ft,  A)  is  given, 
for  k  in  k|  ,k? _ ,k/vk  do 

Find  eigenvalues  and  eigenvectors  B'k  of  H(k.  ft , ft,  ft,  A). 

end  for 


Figure  20.  Pseudocode  for  Luttinger-Kohn  formulation  of  the  k  ■  p 
method. 


4.3  The  Envelope  Function  Approximation 

If  one  substitutes  the  basis  function  expansion  in  equation  109  into  the  Bloch  theorem  expression 
for  the  one-electron  wavefunction  in  equation  10,  one  obtains  the  following  expression  for  the 
one-electron  wavefunction: 

^ik(r)  =  ^[c£k)eikr](/m( r)  =  ^  tfk)(r)(/m(r)  .  (163) 

m  m 

Here,  ?/>ik(r)  is  a  sum  of  modulated  periodic  functions.  The  functions  Um  (r)  are  periodic 
functions  with  the  same  periodicity;  that  is,  Um(r  +  R)  =  Um( r)  where  R  is  a  lattice  vector  of 
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the  crystal.  The  functions  f^k>  (r)  that  modulate  the  amplitudes  of  the  functions  Um(  r)  are 
envelope  functions.  In  general,  they  do  not  share  the  periodicity  of  the  Um  (r),  i.e.,  f^k>  (r) 

A  fmk> (r  +  R),  and  since  k  is  restricted  to  the  first  Brillouin  zone,  f^k> (r)  =  C1^k')elk  r  is  a 
long- wavelength  function  that  varies  far  more  slowly  than  Um  (r).  The  expansion  for  the  one- 
electron  wavefunction  relies  on  Bloch  theorem  and  thus  only  applies  to  systems  that  are  periodic 
in  all  three  spatial  dimensions.  The  general  idea  behind  the  envelope  function  approximation 
(86,  87),  though,  is  that  even  for  systems  that  lack  periodicity  along  one  or  more  directions,  the 
one-electron  wavefunction  may  still  be  approximately  expressed  as  a  sum  of  modulated  periodic 
functions;  that  is, 

Va  (r)  =  ^  (r)  Um'(  r) ,  (]  64 

m' 


where  a  is  a  general  label  for  the  function.  For  a  system  that  has  two-dimensional  periodicity, 
such  as  a  system  composed  of  a  layer  of  one  material  sandwiched  between  two  layers  of  a 
different  material,  the  Bloch  theorem  still  holds  for  planes  normal  to  the  layer,  so  for  a  normal 

direction  along  r3,  the  envelope  function  is  r)  —  fmk]]\r)  —  fm  (r3)elk|l'r|1  where  k||  = 
(k1,k2, 0)  and  rn  =  (rl7r2,0)  (86).  For  a  wire  pointing  along  r3, = /J^tfe3')(r)  = 
fm\ri ,r2)elk3r3  (88).  Otherwise,  there  is  no  Bloch  wave  vector  at  all,  and  the  label  a  is  simply 
the  index  i.  The  envelope  function  is  taken  to  be  slowly  varying  in  comparison  to  Um.  With  the 
envelope  function  equation  in  place,  a  matrix  differential  equation  analogous  to  the  matrix 
equation  1 10  for  the  k  ■  p  method  for  bulk  crystals  may  be  obtained: 


Eafi“\r) 


m 


(165) 


where 


E  FA,  tot 


H m'm 


(r,  Wsys) 


y  0,(r)/?'“({R}L). 


(166) 


and 


-  H^mi) + -fk  •  Prn'm  +  sr, 


h2k2 


<erm({R)D  = 


a 


m  m  1  '-’m  m  2me  ’ 

Um' (r)^ie-(r,  {R}l)um(r)d3r . 


(167) 

(168) 


Here,  k  =  ^  and  k2  —  k  •  k.  For  a  direction  i  along  which  the  system  is  periodic,  k-if^ 

=  kifm  J  where  kL  is  an  element  of  the  Bloch  wavevector.  {R}sys  is  shorthand  for  the  positions 
and  species  of  the  atoms  in  the  system  (which  is,  in  general,  not  equal  to  {R}Lp  the  positions, 
and  species  of  the  atoms  in  a  bulk  crystal  of  type  l.  ®i(r)  is  a  step  function  that  is  1  if  point  r  is 
within  a  material  of  type  l  and  is  zero  otherwise.  Away  from  an  interface,  for  a  point  r  within  a 
material  of  type  Z,  H^tot(r,  {R}sys)  =  ({ R}L>). 
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The  approximations  involved  in  equations  165-167  may  be  made  clear  through  a  brief  outline  of 
their  derivation  by  Burt  (57).  If  one  substitutes  the  envelope  function  expansion,  equation  164, 
into  the  one-electron  Schrodinger  equation,  one  may  obtain  the  following  equation: 

h2k2 


I 


P  P 

2m0 


Um'(X) 


+  ^-[pf/m'(r)]  -k  +  I/^Cr) 


2  mP 


i 

+  ^ext.eff(r,{R}sys)^'(r)  -  Ea U m< (r)j  f%? (r)  =  0  . 


(169) 


Since  derivatives  of  Um( r)  have  the  same  periodicity  as  Um( r)  itself,  one  may  write,  without 
approximation  (57,  57), 


^  **Um'  (r)=X  /  , 


2  mp 


m  c 


2mc 


and 


P  Um'ix) 


(170) 


(171) 


where  the  sums  may  be  infinite.  If  the  effective  potential  can  be  approximated  as 

^ext,eff(r<  Wsys)  ~  ^  Qf  (r)l/ext,cff(r»  WL)  , 


(172) 


then  (57,  57) 


^ext,eff(l7  {R}sys)^a  (0  ~ 

SZfl  ^'(r)Vext,eff(r.  {R}L)l/m(r)d3r 


m  m 


fiat\ r)  ‘ 

Jm  K  J 


(173) 


Once  one  substitutes  equations  170,  171,  and  173  into  equation  169,  the  matrix  differential 
equation  165  readily  follows. 

The  approximation  of  effective  potential  Text, eft  (r,  (Rlsys)  in  equation  172  assumes  the  potential 
field  is  due  to  atoms  in  a  regular  lattice  array.  This  is  reasonable  for  a  point  r  in  an  environment 
that  locally  resembles  a  bulk  crystal,  such  as,  for  example,  in  the  interior  of  a  layer  of  a 
superlattice  composed  of  alternating  single-crystal  layers,  provided  that  the  layers  are  sufficiently 
thick.  However,  for  a  point  r  near  an  interface,  this  approximation  is  rather  crude.  In  actuality, 
the  interface  itself  can  be  considered  a  region  of  small  but  finite  thickness,  as  illustrated  in  figure 
21,  which  shows  an  interface  between  GaAs  and  AlAs.  The  yellow  lines  mark  planes  of  A1  and 
Ga  atoms  on  either  side  of  the  interface,  with  a  plane  of  As  atoms  in  between.  Across  the  region 
between  the  planes  of  Ga  and  A1  atoms,  the  true  potential  felt  by  the  electron  (i.e.,  equations  5 
and  6)  is  generally  continuous,  though  it  may  rise  or  fall  steeply  across  this  interface. 
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Figure  21.  Interface  between  AlAs  at  top  (grey 
and  purple  atoms)  and  GaAs  at 
bottom  (pink  and  purple  atoms). 

The  two  yellow  dotted  lines  indicate 
the  planes  of  A1  and  Ga  atoms  on 
opposite  sides  of  the  interface.  The 
figure  was  created  with  Jmol  (45). 

Since  the  envelope  function  expansion  contains,  in  principle,  an  infinite  number  of  terms,  the 
matrix  differential  operator  (R}sys)  likewise  has  an  infinite  number  of  rows  and 

columns.  One  way  of  obtaining  a  finite  matrix  operator  from  is  to  note  how 

^m'm ((R}oo)  in  equation  167  and  //^^({R}^)  in  equation  111  are  analogous,  with  the  former 
obtainable  from  the  latter  by  making  the  replacement  k  ->  k.  This  replacement,  then,  may  be 
made  for  a  previously  fonnulated  k  ■  p  Hamiltonian  matrix  for  a  bulk  crystal  (86).  However, 
there  are  problems  with  this  approach.  Material  parameters  are  functions  of  position  across  an 
interface,  so,  for  example,  the  Luttinger  parameters  would  become  ym( r)  =  X;  ®i(x)Ym  where 
Ym  is  the  mth  Luttinger  parameter  for  material  l.  While  the  scalar  components  of  the  Bloch 
wave  vector  k  commute  with  such  parameters,  components  of  the  differential  operator  k  do  not 

(89) .  A  work-around  for  this  issue  has  been  to  impose  a  symmetrization  scheme  on  /7^^tot 

(90) ,  so  that,  for  example,  operators  within  it  of  the  forms 

q2  ^ 

4(x)— andS(x)-,  (174) 

ox1  ox 

are  replaced,  respectively,  by 

d  d  if  d  d  \ 

a„d_^W_  +  _eWj.  (175) 

This  operator  replacement  scheme,  however,  is  ad  hoc  (91)  and  can  lead  to  spurious  solutions 
(92). 

Instead  of  directly  substituting  k  for  k  in  a  previously  formulated  finite-sized  k  ■  p  Hamiltonian 
matrix,  matrix  differential  equation  165  can  itself  be  transformed  into  a  matrix  equation  with  a 
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finite  number  of  terms,  following  the  method  of  Burt  (87).  In  this  method,  the  expansion  in 
equation  164  is  split  into  two  sums, 

■Mr)  =  Zifs(a\r)UM  +  E? f™(r)l/r(r) .  (176) 

Following  the  notation  given  previously  for  the  derivation  of  the  Luttinger-Kohn  Hamiltonian 
matrix,  the  states  of  the  first  sum  are  in  class  A,  and  the  states  in  the  second  sum  are  in  class  B . 
The  class  A  states  are  not  necessarily  the  exact  same  ones  used  by  Luttinger  and  Kohn,  but  they 
are  the  states  that  contribute  the  most  to  xpa(r).  The  index  s  refers  to  envelope  functions  and 
tenns  corresponding  to  class  A  states;  index  r  refers  to  envelope  functions  and  tenns 
corresponding  to  class  B  states,  and  index  m  refers  to  terms  in  either  sum.  With  this  notation, 
the  envelope  function  matrix  equations  may  be  expressed  as 


V  _l_  JL u  .  n  4-  X  h2k2)  _  p  f(a)rr^ 

/-im  [nrmyl  J  Vrm  '  urm  2jn  Jm  v*  J  Lla)r  V*  J  ■> 


(177) 


and 


where 


y  fi(r)  +—  kpsrn  +  5: 
Z_i  mP 


h2k 2 


2  mP 


fm\ r)  =  £»//“’( r) , 


(178) 


(179) 


Burt  neglects  the  curvature  of  the  class  B  envelope  functions  k2f^\ r)  is,  so  equation  177 
becomes  approximately  (57) 


Z  +  Z  K' «  +  J-k  •  P 

r'  i,s'  6 


r)  *  S»/"'(r)  , 


(a). 


(180) 


where  the  sum  over  m  has  been  split  into  two  sums  over  r'  and  s',  where  r'  and  s'  correspond  to 
classes  A  and  B,  and  the  terms  h  k  •  p  rr'  /me  are  also  assumed  to  be  negligibly  small  (57).  If  the 
class  B  periodic  functions  I/r(r)  are  such  that  77^’J  =  H^?8rrr,  then  the  class  B  envelope 
functions  are  approximately  (57,  57) 

1 


/r  a)  (r) 


Ea-H%\ r)V 


Yte(r)+-  k-p 

Z_i  L  rs 


(181) 


Equation  178  can  be  rewritten  to  show  its  dependence  on  /r  (r): 
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(182) 


r 


Equation  181  is  then  substituted  into  equation  182,  but  several  of  the  resulting  terms  are  either 
zero  due  to  symmetries  in  class  A  states  or  are  neglected  because  they  are  only  significant  near 
an  interface.  (Terms  containing  derivatives  of  matrix  element  Hrs  are  an  example  of  the  latter, 
since  these  elements  are  constant  within  a  bulk  material,  and  regions  away  from  the  interface  are 
presumed  to  be  bulk-like.)  This  yields  the  following  (81,  87): 


=  £0//“V) 


r(a)/ 
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where  r,  {R)sys)  is  called  a  Burt-Foreman  Hamiltonian.  Here,  psr  is  position-independent 
and  commutes  with  k.  (r),  however,  is  clearly  not  position- independent.  A  Hamiltonian 
matrix  derived  from  equations  183-185  should  have  the  proper  ordering  between  the  operator  k 
and  the  empirical  parameters  formed  from  the  sums  over  r  (class  B).  One  can  see  an  example  of 
such  a  Hamiltonian  from  Foreman  (93).  If  Ea  ~  E0  and  =  Em08mim,  then  for  a  bulk 

material,  H^i  resembles  the  Luttinger-Kohn  Hamiltonian  shown  in  equation  157. 


Figure  22  is  pseudocode  that  shows  the  outline  of  how  an  implementation  of  the  envelope 
function  approximation  may  be  done.  In  this  pseudocode,  discretization  refers  to  the  use  of 
either  finite  difference  or  finite-element  methods  to  transform  a  matrix  differential  eigenequation 
Hf  A)(r)  =  Eia f(a)(r),  where  the  elements  of  H  are  differential  operators  and  the  elements  of 
fa(r)  are  functions,  into  a  matrix  eigenvalue  problem,  HF1'"-’  =  Eia¥^a\  where  the  elements  of 
H  and  F (ai  are  numerical  values. 
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Set  of  empirical  parameters  params  is  given. 

Matrix  differential  operator  H(k\,  fa,  fo, params)  is  given. 

Number  of  dimensions  along  which  electrons  are  confined.  No,  is  given. 

if  No  =  I  then  {Electrons  are  free  to  move  within  a  plane.} 

Let  differential  operators  k\  and  £2  become  the  scalars  k\  and  ki. 

Discretize  the  matrix  differential  equation  H(£|  ,k2,k3,params)tf'k]\r)  =  £,*  f^'k  *(r), 
where  k  =  (k\,k2,0),  along  a  line  pointing  along  the  axis  (i.e.  ^-direction),  forming  a 
matrix  eigenvalue  problem  H(k  =  £,k  Fl'kl^. 

for  k\  from  to  k™'dX  do 
for  ki  from  kfm  to  /c™*  do 
Let  k||  =  (k\,k2,0). 

Find  eigenvalues  £,* t  and  eigenvectors  F!,k  1  of  H(k|| ). 

end  for 
end  for 

else  if  No  =  2  then  {Electrons  can  only  move  along  line.} 

Let  differential  operator  k 3  become  the  scalar  k$. 

Discretize  the  matrix  differential  equation  H(&i  ,k2,ks, params) f^'*3^ (r)  =  £’,Jt,f^'*3^(r) 

in  the  plane  (i.e.  xy  plane),  forming  a  matrix  eigenvalue  problem  H ( kj ) F( 1/13 *  =  £,r.  E1'*-'*. 

for  k3  from  kfn  to  k"ux  do 

Find  eigenvalues  E^ ,  and  eigenvectors  FI'*3)  of  H^). 

end  for 

else  if  No  =  3  then  {Electrons  confined  in  all  three  dimensions.} 

Discretize  the  matrix  differential  equation  H(k\ , fa, params)  (r)  =  £,f^(r)  in 
all  three  dimensions,  forming  a  matrix  HF*'*  =  £,F''). 

Find  eigenvalues  E,  and  eigenvectors  F(,)  of  H. 

end  if 

Figure  22.  Pseudocode  for  implementing  the  envelope  function  approximation.  The  matrix  differential 

operator  H  (k  {lk  ,k3,  params)  may  either  be  a  symmetrized  bulk  k  ■  p  Hamiltonian  matrix 
where  the  substitution  k  ->  k  has  been  made  or  a  Burt-Foreman  Hamiltonian. 
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The  envelope  function  expansion  of  the  one-electron  wavefunction  in  equation  164  is  not  fully 
general.  In  a  region  of  a  structure  that  is  locally  periodic — that  is,  a  region  composed  of  unit 
cells  that  may  be  tiled  throughout  a  finite  space  rather  than  to  infinity — it  is  physically 
reasonable.  The  local  periodicity  is  captured  by  the  periodic  functions  Um  (r) ,  while  the  lack  of 
full  periodicity  due  to  the  presence  of  features  such  as  material  interfaces  in  the  structure  is 
accounted  for  by  the  envelope  functions  .  However,  regions  of  structures  may  lack  even 
approximate  local  periodicity,  for  example,  if  they  are  sufficiently  inhomogeneously  strained. 
(Such  regions  may  also  be  regions  where  the  potential  cannot  be  approximated  as  bulk-like.) 

The  significance  of  this  can  be  seen  in  figure  23,  where  Um(  r)  is  superimposed  over  an 
undistorted  lattice  with  the  same  periodicity  as  Um( r)  and  a  distorted  lattice  that  lacks 
periodicity.  A  slowly  varying  envelope  function  can  modify  the  amplitude  of  Um(r),  but  it 
cannot  stretch  or  compress  parts  of  it  to  match  the  distortions  in  the  lattice.  Also,  unlike  the 
periodic  functions  in  Fourier  series  expansions,  the  periodic  functions  in  the  envelope  function 
expansion  of  equation  164  all  have  the  same  periodicity  over  the  whole  domain  CL.  Lattice 
mismatch  in  a  structure,  though,  means  that  different  regions  within  the  structure  may  have 
different  local  periodicities,  so  this  would  limit  the  envelope  function  expansion — at  least  the  one 
given  in  equation  164 — to  lattice-matched  structures,  such  as  those  made  of  GaAs  and  AlAs. 


Figure  23.  Schematic  of  a  periodic 
function  Un( r) 
superimposed  over  a  lattice 
of  evenly  spaced  atoms  (top) 
and  a  distorted  lattice  where 
the  spacing  between 
neighboring  atoms  is  no 
longer  the  same. 

However,  the  envelope  function  approximation  has,  nonetheless,  been  applied  to  structures  with 
significant  lattice  mismatch,  such  as  a  gallium  nitride  (GaN)  dot  embedded  in  an  AIN  matrix 
(94),  or  an  InAs  quantum  dot  embedded  in  a  GaAs  matrix  (89).  Work  by  Foreman  (95)  indicates 
why  this  appears  to  lead  to  reasonable  results  in  practice.  They  show  the  one-electron 
wavefunction  may  be  expanded  as 


</>«(  r) 


Zm'  If-""  ?“( Ri)aB(r  -  R«)Um'(r,  Rj)  , 


(186) 
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where  Nce\\  is  the  number  of  crystal  unit  cells  in  the  system  (which  in  general  may  not  be 
identical  when  the  system  is  composed  of  multiple  materials),  R,  is  the  position  vector  pointing 
to  unit  cell  i,  and  SB(r  —  Rt)  is  defined  as 


SB(r-Ri) 


£li 

(27t)3  jFBZ 


f  eikrd3r 

JFR7  c  u  1  ■ 


(187) 


Here,  “FBZ”  is  the  first  Brillouin  zone  of  a  crystal  unit  cell  at  Rj,  and  H,  is  the  volume  of  this 
cell.  This  function  is  similar  to  the  Dirac  delta  function  in  that  it  peaks  sharply  near  Rj,  and  is 
zero  outside  a  neighborhood  of  Rj,  with  a  volume  about  equal  to  flj.  Um<  (r,  Rt)  is  the  value  of 
Um'(r )  for  a  bulk  crystal  whose  unit  cell  is  the  same  as  that  of  cell  i.  F^/(Ri)  is  a  discrete 
envelope  function  defined  at  points  Rj,  and  its  relationship  to  the  continuous  envelope  function  is 

C\ r)  =  Sfe"  -  R,)  •  (188) 

If  Umt  (r,  Rj)  is  the  same  for  all  crystal  unit  cells  in  the  system,  then  Umt (r,  Rj)  =  Umt  (r),  and 
the  expansion  in  equation  186  reduces  to  that  in  equation  164.  However,  this  expansion  is  more 
general  than  the  original  envelope  function  expansion  and  accounts  for  changes  in  the  local 
periodicity.  With  this  redefinition  of  fff( r),  the  matrix  operator  (r,  (R}sys)  in  equation 

165  becomes 
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<t(r)  =  S^Kr)^;({R}L) 


(190) 


Pm'mW  =  Si0i(r)/n.H^,(r,Rj)pI/rn(r,Rj)d3r,  (191) 

pm'm(r)  =Si0i(r)/n.I/^,(r;Rj)p[5iJ(r-Rj)t/7n(r;Rj)]d3r,  (192) 
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(193) 


1,  r  in  unit  cell  i 
0,  otherwise 


Here,  {R}L  is  the  set  of  positions  and  species  of  the  atoms  of  a  bulk  crystal  whose  unit  cell  is  the 
same  as  that  of  cell  i.  The  symmetry  of  the  sum  p^p^Cr)  •  k  +  k  •  p^p^Cr)  in  the  second  term 
of  equation  189  is  not  due  to  an  artificially  imposed  symmetrization  scheme,  such  as  the  one 
discussed  above  and  shown  in  equation  175.  If  the  structure  in  question  is  composed  of  locally 
periodic  crystal  regions  separated  by  material  interfaces,  then  away  from  an  interface,  p m>m  is 
effectively  independent  of  r,  Pm/m(r)  is  zero,  and  H^,^ot(r,  (R}sys)  in  equation  189  reduces  to 
the  /?“?"0,( r,  {R}sys)  in  equation  167.  The  tenns  containing  Pm/m(r),  then,  are  interfacial 
terms,  which  are  often  negligible,  especially  in  systems  composed  predominantly  of  regions  of 
bulk  crystal,  e.g.,  superlattices  with  thick  layers  (57).  The  end  result  is  that  the  envelope 
function  matrix  equation,  i.e.,  equation  165,  may  be  applied  to  structures  that  are  not  necessarily 
lattice-matched. 


4.4  Effective  Mass  Approximation 


Sometimes  the  tenn  “effective  mass  approximation”  is  used  as  a  synonym  for  the  k  ■  p  and 
envelope  function  methods  (32,  82,  87),  but  here  it  will  be  used  to  describe  a  particular 
simplification  of  these  methods.  One  can  see  this  approximation  in  the  solutions  to  the  Kane 
Hamiltonian  matrix  for  small  k,  equations  141,  which  all  have  the  following  form: 


Eik  -  Ei  + 


h2k2 
2  m; 


(194) 


This  is  a  solution  to  the  Schrodinger  equation  for  an  independent  particle  with  mass  mL  in  a 
potential  with  the  constant  value  EL, 

+  Ei\  =  Eikipik(r)  ,  (195) 


and  because  of  this,  is  called  the  effective  mass.  Accordingly,  electrons  with  a  small  Bloch 

wave  vector  k  behave  approximately  as  if  they  were  free  electrons  with  this  effective  mass. 
Holes  generally  behave  approximately  as  if  they  were  free  electrons  with  a  negative  effective 
mass. 


In  the  envelope  function  method,  the  effective  mass  approximation  amounts  to  an  envelope 
function  expansion  with  only  one  term,  that  is, 

■Mr)  =  /(«>(r)(/(r)  ,  (196) 

and  the  envelope  function  matrix  differential  equation  reduces  to  the  following  one-particle 
equation  (10, 17), 
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Strictly  speaking,  equation  197  is  valid  only  if  mi  is  constant.  In  a  heterostructure,  is  at  most 
piecewise  constant,  taking  on  one  value  in  a  part  of  the  structure  composed  of  one  material  and  a 
different  value  in  a  part  composed  of  a  different  material.  (Similarly,  £)  can  at  most  be  only 
piecewise  constant.)  In  such  a  structure,  then,  mi  varies  with  r,  so  the  operator  p  •  p/2m(  is  no 
longer  Hermitian,  so  Ea  is  no  longer  guaranteed  to  be  real.  This  problem  can  be  avoided  by 
modifying  the  equation  to  become  the  following  (10): 

[p  ‘  (l^yj)  +  E] /(a)(0  =  EJ^(r)  ,  (198) 

or  equivalently, 

[k  •  (^)  +  £,] /<“>( r)  =  £«/w( r)  .  (199) 

Figure  24  shows  a  pseudocode  outline  of  an  implementation  of  the  effective  mass  approximation 
for  the  case  where  the  effective  mass  is  either  piecewise  constant  or,  for  the  pure  bulk  case, 
constant.  In  this  pseudocode,  0;(r)  is  a  step  function  that  is  1  if  point  r  is  within  a  material  of 
type  Z  and  is  zero  otherwise.  Here,  discretization  refers  to  the  use  of  either  finite  difference  or 
finite  element  methods  to  transfonn  a  differential  eigenequation  //f (r)  =  Eaf^(r)  into  a 
numerical  matrix  eigenvalue  problem  HF("Z  =  Ea  F (ai . 

One  use  of  the  effective  mass  approximation  is  to  supplement  a  k  ■  p  or  envelope  function 
implementation  that  only  solves  for  the  energies  of  valence  bands.  That  is,  if  one  uses  a 
Luttinger-Kohn  Hamiltonian  to  determine  the  heavy-  and  light-hole  bands,  one  may  then  use  the 
effective  mass  approximation  to  estimate  the  behavior  of  the  conduction  band  near  k  =  0. 

4.5  Accounting  for  Strain  in  the  k  •  p  and  Envelope  Function  Methods 

In  ab  initio  and  empirical  atomistic  methods,  the  effects  of  strain  are  at  least  partially  taken  into 
account  through  the  atomic  positions  (R)M  that  are  the  input  of  the  effective  potential 
Fext,cff(r>  (R)„).  However,  in  the  k  ■  p  method,  the  effect  of  I4xteff  is  taken  into  account 
indirectly  through  parameters  such  as  Eg,  A,  the  Kane  parameter  PK ane,  or  the  Luttinger 
parameters.  Atomic  positions  are  not  an  input  to  the  method  at  all.  Instead,  a  strain-dependent 
tenn  is  added  to  the  effective  potential  operator,  usually  using  the  approach  of  Bir  and  Pikus  ( 1 ). 


65 


For  each  material  /,  effective  mass  m,(/)  is  given. 

For  each  material  /,  band  extremum  £,(/)  is  given. 

Number  of  dimensions  along  which  electrons  are  confined,  Np,  is  given. 

if  Np  =  0  then  {Electrons  are  free  to  move  in  any  direction.} 

Let  m,  =  m,(l)  and  £,  =  £,(/)  (There  is  only  one  bulk  material  /  here.} 

Energy  band  as  a  function  of  wavevector  k  is  £,k  Ej  +  /i2|k|2/2m,\ 

else 

Let  Hamiltonian  operator  H(ic\, ki, kj)  «  J^0/(r)  (  k  -  —jjr  +  £,(/)),  where  k  =  (ki.ic2,k},). 

if  Np  =  I  then  {Electrons  are  free  to  move  within  a  plane.} 

Let  differential  operators  £|  and  £2  become  the  scalars  k \  and  Zi- 

Discretize  the  differential  equation  fi{k\,ki,k-i) /(,k|l*(r)  =  £,k|/*‘k|  ’(r), 
where  k  =  ( k\,k2,0 ),  along  a  line  pointing  along  the  r-$  axis  (i.e.  --direction),  forming  a 
matrix  eigenvalue  problem  H(k||)F*,kH*  =  £,k 

for  k  1  from  k""n  to  k"VM  do 
for  k2  from  kfn  to  do 
Let  k  =  ( k\  ,k2,0). 

Find  eigenvalues  £,k  and  eigenvectors  Fl,kn 1  of  H(k  ). 

end  for 
end  for 

else  if  Np  =  2  then  {Electrons  can  only  move  along  line.} 

Let  differential  operator  £3  become  the  scalar  £3. 

Discretize  the  differential  equation  //(Z|  /*'*’*(r)  =  £,*, /*'*’'( r) 

in  the  rir?  plane  (i.e.  x y  plane),  forming  a  matrix  eigenvalue  problem  H(A'3)F!,<,,)  =  E/^F*'^. 

for  ky  from  kfn  to  k"VM  do 

Find  eigenvalues  £,•*,  and  eigenvectors  F*'*3*  of  H (Z:^ ) . 

end  for 

else  if  Np  =  3  then  {Electrons  confined  in  all  three  dimensions.} 

Discretize  the  differential  equation  H(k\,k2,ki)  =  £,/*'' (r)  in 

all  three  dimensions,  forming  a  matrix  HF(,)  =  £,F^. 

Find  eigenvalues  Ej  and  eigenvectors  F1'1  of  H. 

end  if 
end  if 


Figure  24.  Pseudocode  for  outline  of  implementation  of  effective  mass  approximation  for  (piecewise) 
constant  effective  mass. 


66 


The  positions  and  species  of  atoms  in  a  strained  bulk  crystal  will  be  denoted  here  as  {R}^,  and, 
accordingly,  the  one-electron  Schrodinger  equation  for  a  strained  bulk  crystal  where  spin-orbit 
coupling  is  ignored  is 

Hie- (r,  W*  )</hk (r)  =  [^  +  kext,eff(r.  {R}*  )]  ipt k(r)  =  E?kipik( r) ,  (200) 

where  Efk  denotes  an  eigenvalue  for  the  strained  system.  The  one-electron  Hamiltonian  for  the 
strained  crystal  may  be  rewritten  as 

#ie-(r,{R}*)  =  tfle-(r,{  R}oo)  +  [kext,eff(r,{R}^)  -  tUeff(r.{R}oo)] ,  (201) 

where  (R}oo  here  denotes  the  positions  and  species  of  atoms  in  an  unstrained  bulk  crystal.  The 
term  kext,eff(r,  ml)  —  kext,eff(r»  (Rjoo),  though,  cannot  be  regarded  as  a  small  perturbation  on 
the  Hamiltonian  for  the  unstrained  system,  Hle-  (r,  (R)oo),  because,  as  Bir  and  Pikus  point  out 
(/),  this  tenn  is  not  generally  small.  This  can  be  illustrated  in  figure  25,  which  shows  model 
undeformed  and  deformed  sinusoidal  potentials  of  a  Active  one-dimensional  crystal,  denoted  v 
and  vx,  respectively,  and  the  differences  between  them.  A  small  uniform  strain  e  is  introduced 
that  stretches  the  period  of  v  by  a  factor  of  1  +  e  and  changes  its  amplitude  slightly  by  we, 
where  w  is  a  parameter.  Although  these  model  potentials  have  only  slightly  different  periods, 
with  v(x )  =  v0  sin  x  and  vx(x)  —  ( v0  +  we)sin[x/(l  +  e)],  the  difference  between  them 
grows  as  x  increases  to  become  on  the  same  order  as  the  potentials  themselves. 


Figure  25.  Undeformed  and  deformed  sinusoidal  potentials  of  a  one-dimensional  crystal, 
and  the  difference  between  them.  v(x)  =  v0  sin  x,  and  v  (x) 

=  (v0  +  we)  sin  [x/(l  +  e)].  For  the  graph  above,  v0  =  1,  e  =  1%,  and  w  =  -0.3. 
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In  the  previous  one-dimensional  example,  one  can  specify  a  deformation  map  x  that  relates  the 
coordinate  of  the  unstrained  system,  now  denoted  x°,  to  the  coordinate  of  the  strained  system; 
that  is,  x  =  x(x° )  =  (1  +  e)x° .  One  can  also  observe  that  while  the  difference  vx(x)  —  v(x)  is 
not  small,  vx(x(x°))  ~  v(x°)  —  we  is  small,  provided  that  the  parameter  w  is  not  too  large. 
Similarly,  the  difference 

Vext,eff(x(r°),  (R&)  -  Vext,eff(r°,  Woo)  *  ^(r°- Woo >ij  ,  (202) 

where  r°  is  the  coordinate  vector  for  the  unstrained  system,  is  also  small,  provided  that  VLj  is  not 
too  large.  Again,  the  convention  of  summing  over  repeated  indices  is  used.  For  a  bulk  crystal 
subjected  to  homogeneous  strain,  the  deformation  map  x(r°)  equals  Fr°  =  (I  +  V0u)r°,  where 
F  is  the  deformation  gradient  and  V0u  is  the  gradient  of  the  displacement  with  respect  to  the 
coordinates  of  the  unstrained  system.  If  there  is  no  rotation,  the  displacement  gradient  is 
symmetric  and  the  small-strain  tensor  may  be  taken  to  be  e  =  V0u  (66),  and  x(r°)  =  (I  +  e)r°. 
(Use  of  the  small  strain  tensor,  however,  implies  that  V0u  «  Vu  where  Vu  is  the  gradient  of  the 
displacement  with  respect  to  the  coordinates  of  the  strained  system  [96].)  The  coordinate 
transformation  r  =  /(r°)  is  now  introduced  into  the  one-electron  Schrodinger  equation  for  a 
strained  crystal: 

Hie- (x(r°),  W*  )</hk(x(r0))  =  E?kxp  ik(x(r®))  .  (203) 

The  momentum  operator  within  Hle~  (x(r°)<  {R}^)  needs  to  be  expressed  in  tenns  of  the 
coordinates  of  the  unstrained  system  (32).  Applying  the  chain  rule  yields 


Pi  = 


dr®  d 


dr® 


(-ift) 


dr® 


Pi 


(204) 


where  p°  =  (p®,  p®,  p®)  is  the  momentum  operator  of  the  one-electron  Hamiltonian  of  the 
unstrained  system  in  terms  of  r°.  For  small  strain,  the  inverse  deformation  map  r°  =  X-1(r)  is 
approximately  (32, 1) 


X-1(r)  =  (I-€)r,  (205) 

and,  accordingly, 

Pi  -  fa  -  e^p®  .  (206) 


Once  the  momentum  operator  is  expressed  in  tenns  of  unstrained  coordinates,  then  one  may 
write 
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P  P  =  Pipi  ~  (Sij  -  £ij)pj(8ik  -  eik)v°k  , 

—  {8jk  —  2  6 jk  +  £ij£jk)PjPk  > 
and 

-  -  2 €jk)pyk  ,  (207) 

where  the  small  strain  assumption  means  that  the  product  of  strain  components  is  negligible. 

The  one-electron  Hamiltonian,  then,  becomes  (32, 7 ) 


ffie-(x(r°),{R}*)  =  Hle-( r°,{R}oo)  +  eJkSJk  , 

where 

_  p°  •  p° 

Hle-  (r°,  (Rjoo)  =  +  Kext,eff(r°,  {R}^)  , 

2m  e 

and 

rue 


(208) 


(209) 

(210) 


The  one-electron  wavefunction  satisfies  the  Bloch  theorem  in  both  strained  and  unstrained 
coordinates.  To  show  this,  the  Bloch  wave  vector  k  needs  to  be  expressed  in  terms  of  unstrained 
coordinates.  The  primitive  lattice  vectors  and  reciprocal  lattice  vectors  in  both  strained  and 
unstrained  coordinates  still  satisfy  the  following  relationships  (72), 

•  a j  —  2n8ij,  and  b°  •  a °  =  2n8ij  ,  (211) 

where  the  superscript  “0”  indicates  that  a  quantity  pertains  to  the  unstrained  system.  Since 
a j  —  (I  +  e)a j* ,  b*  =  (I  +  e)-1b°  ,  and 


k  =  ^  b^  =  (I  +  e)_1  ^  b°  =  (I  +  e)_1k° 


(212) 


i=i 


i=i 

Jkr  _  „ik°r° 


where  is  a  real  number.  Accordingly,  elk  r  =  elk  r  ,  and 

^ik(r)  =  eikruik(r)  =  eik°'r°ui(I+e)-iko((I  +  e)r°) 
-eik°rVk0(r°)=^fk0(r0), 

where  u^k0(r°)  =  Ui(I+e)-iko((I  +  e)r°).  The  function  ufko(r°)  has  the  periodicity  of  the 
unstrained  lattice  (7);  that  is, 


(213) 
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and 


Ufk°(r°  +  =  Ut(I+e)-1k°((I  +  €)(r°  +  R°))  5 

—  utk((I  +  e)(r°  +  R°))  5 

—  utk((I  +  e)r°  +  (I  +  e)R°)  5 
=  wik(r  +  R)  =  Mik(r)  , 


=  wi(I+e)-1k°((I  +  £)r°)  =  Ufk°  (r°)’ 

where  R  and  R°  are  lattice  vectors  of  the  strained  and  unstrained  crystal. 
Applying  the  operator  £jkDjk  to  ip*ko(r°)  yields 


ejkDjk4>nP  (r°)  = 


ik°-r° 

[tjkPjPk  +  CjkihtfpZ  +  hk°kp°k)  +  Ejkh2  kk\u*ko  (r°) 
+  f;7c^7c(r0,{R)oo)wfk0^r°)  ’ 


and 


ik°-r° 

--ZI—l€jkPjPk  +  2ejkhkfp°k  +  ejkh2k°k%\u*k0(r0) 

Trie 

+  ^7c^7c(r0,{R}oo)wfko(r°) > 


ik°-r° 

=  elk  r  e 


Jk 


Jjk 


2hkjpk 

o^; 

o  — 

CN 

■*; 

1 

me 

me  J 

Ufk°(r°)  1 


This  leads  to  the  following  Schrodinger-like  equation: 

Hie-  (r°,{R}00)  +  ^-{8Jk  ~  2 ejk)kjpk  +  ejk Djk  ufk0(: r°) 

Ule  J 

Eik°-^-  i5jk  -  2ejk)kfk%\  wfko(r°)  . 


(214) 


(215) 


(216) 


As  before,  the  periodic  part  of  the  wavefunction  is  expanded  in  terms  of  the  states  where  the 
Bloch  wave  vector  is  zero,  i.e., 


wfko(r°)  =  ^  ci'k  ]Um(r) 

m 

and  a  matrix  equation  is  obtained: 


h2 


(«,*  -  2elk)k?ki 


/Aik) 
Lm'  ’ 


m  L 

R)- <0  =  AVMJ  +  £ fe  -  2eit)kfp°X  + 
Pm'm  =  f 


O'k) 

m'm 


(217) 


(218) 

(219) 

(220) 
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and 


Dm')m=[  ul'WDijUmi r)d3r,  (221) 

Jnc 

with  ac  being  the  volume  of  the  unit  cell  of  the  crystal.  Since  the  strain  here  is  taken  to  be 
homogeneous,  e is  a  constant  and  can  be  moved  outside  of  any  integration. 

When  the  basis  functions  Um  (r)  =  um0  (r)  are  the  Luttinger-Kohn  basis  functions  in  equation 
145,  then  ejkD(jl,c>  is  such  that  (32) 

J  ^  771  771  v  7 


D 


O'fc) " 

i'm 
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Se 
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-st 

Pe-Qe 

0 

Pt 

0 

Pe-Qe 
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Pt 

s t 
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4V2  Se 

V2  R\ 

4V2St 

V2  Qt 

0 

—  Se/y[2 
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Re 

-V2  Qe 

S£ 

Se 

JV2S'l 

V2  Qe 

Pe  +  Qe 

-V2  Rl 

-St/yfl 

-V2  Re 

Pe 

0 

-S£/V 2 

0 

Pe  - 

where 


and 


P£  -  av  Tr  e  , 
b 

Qe  —  2  (en  +  e22  ~  2e33)  , 


V3 


Re  =  --^On  -e22)  +  ide12. 


(222) 


Se  —  d(e13  —  ie23)  ■  (223) 

The  empirical  parameters  av,  b,  and  d  are  linear  combinations  of  the  matrix  elements  ,  and 
they  related  to  certain  strain  states  (81).  If  the  strain  is  purely  hydrostatic,  then  the  valence  band 
maximum  is  shifted  by 


A  Ev  =  av  Tr  e  . 


(224) 


If  the  strain  is  biaxial  with  en  =  e22  ^  £33,  the  maxima  of  the  heavy  and  light  hole  valence  band 
are  no  longer  the  same,  but  are  separated  by 

AEhi  =  2\be33\ .  (225) 

If  the  strain  is  pure  shear  with  ey  —  es(l  —  SLj),  the  separation  becomes 

AEhi  =  2\des\.  (226) 
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The  conduction  band,  which  is  not  taken  into  account  by  the  strain  correction  to  the  Luttinger- 
Kohn  Hamiltonian  matrix  in  equation  222,  also  shifts  due  to  strain.  For  the  case  of  hydrostatic 
strain,  this  shift  is 

A Ec  =  ac  Tr  e ,  (227) 


where  (32) 

ac  =  [  S+(r)Z\,S(r)d3r  ,  (228) 

Jnc 

and  S  (r)  has  the  same  definition  that  it  does  in  the  formulation  of  the  Kane  k  ■  p  Hamiltonian 
matrix.  A  diagram  illustrating  the  effects  of  strain  on  a  bulk  crystal  band  structure  is  shown  in 
figure  26.  In  general,  compressive  strains  tend  to  increase  the  size  of  the  band  gap,  tensile  strains 
tend  to  decrease  it,  and  strains  departing  from  a  pure  hydrostatic  state  tend  to  separate  the  heavy 
and  light  hole  bands  (97). 


Figure  26.  Schematic  band  structure  of  a  typical  bulk  semiconductor  with  a 
diamond  or  zincblende  crystal  structure,  where  the  solid  lines 
indicate  the  band  structure  of  a  strained  semiconductor,  while  the 
dotted  lines  indicate  the  original  band  structure  before  the 
strain  is  applied.  The  split-off  band  is  not  shown.  A £/,/  is  the 
difference  between  the  maxima  of  the  heavy  hole  and  light  hole 
bands,  which  are  not  necessarily  the  same  once  the 
semiconductor  is  strained.  A Ec  is  the  shift  of  the  conduction 
band  minimum,  while  A Er  ±  A£),//2  is  the  downward  shift  of  the 
light  and  heavy  hole  bands,  respectively  (81).  The  left  and  right 
halves  of  the  horizontal  axis  indicates  the  magnitude  of  k-values 
pointing  along  certain  crystal  directions.  (In  this  schematic,  the 
actual  directions  are  not  important.) 

At  this  point,  a  strain-dependent  term  has  only  been  provided  for  the  one-electron  Schrodinger 
equation  for  an  infinite  periodic  system.  The  empirical  parameters  ac,  av,  b,  and  d  all  pertain  to 
a  bulk  crystal.  Furthermore,  the  strain  is  taken  to  be  homogeneous,  and  the  coordinate 
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transformation  r  =  (I  +  e)r°  implies  that  the  deformation  has  no  rotational  component.  (See  the 
discussion  in  appendix  B  of  a  similar  approximation  in  the  atomistic  strain  formulation  of  Pryor 
et  al.  [95].)  Nonetheless,  when  the  envelope  function  method  is  applied  to  strained  structures, 
usually  what  is  done  is  to  take  a  preexisting  k  ■  p  Hamiltonian  matrix,  substitute  k°  (i.e.,  p °/h) 
for  k°  (possibly  with  an  additional  correction  to  account  for  Burt-Foreman  operator  ordering), 
and  then  add  the  same  strain  correction  term  that  one  would  use  for  a  bulk  k  ■  p  matrix  (10,  32, 
81,  89). 

Blount  (99)  and  Sham  and  Ziman  (100)  have  shown  how  a  strain-dependent  term  may  be 
introduced  into  the  one-electron  Schrodinger  equation  through  a  coordinate  transfonnation  that 
does  not  imply  homogeneity  of  the  strain.  The  coordinate  transformation  is  simply  the 
relationship  between  the  strained  coordinate,  the  unstrained  coordinate,  and  the  displacement  u, 
i.e., 


r  =  x(r°)  =  r°  +  u . 

(229) 

The  momentum  operator  can  be  expressed  as  (99) 

P'  =  (n  -  w)]  p°  = 

(*«-£)>?. 

(230) 

and  so  for  small  strain  (100) 

P  P  =  PiPi  ~  (Sjk  -2^)  p°p£  -  1 

(231) 

Furthermore, 


7  d22± -o  _  ^0  -o  ,  ^0  -o  _  d_2±~ o-o  ,  ^0 

^  drk  PJ  Pk  ~  drk  P1  Pk  +  drk  P1  Pk  ~  drk  P1  Pk  +  drj  PkPJ 

-  drk  pj  Pk  +  drj  pj  Pk  -  +  dr j)  V1  Pk  ~  Z£JkPl  Pk 


duk  __  „  ... 

gr.  eki  +  " ki  > 

where  ooki  is  the  infinitesimal  rotation  tensor  (96).  Therefore, 

P  P  *  (8Jk  -  2 ejk)p0jp°k  -  [p-(eki  +  a>ki)]p% 


( Sjk  ~  2 ejk)PjPk  “  [^rX£ki  +  <Ufci)]  (-ift)pje' 


(232) 

(233) 

(234) 


If  both  eki  and  a)ki  are  sufficiently  slowly  varying,  then  the  second  tenn  in  the  above  equation  is 
negligible,  and  the  remainder  of  the  above  equation  coincides  with  the  corresponding  result  from 
Bir  and  Pikus.  Also,  equations  232  through  234  depend  on  the  assumption  that  the  small  strain 
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tensor  is  formed  from  the  displacement  gradient  with  respect  to  the  coordinates  of  the  strained 
system,  duk/dri  or  Vu,  while  the  fonnulation  of  Bir  and  Pikus  is  derived  under  the  assumption 
that  the  small  strain  tensor  is  fonned  from  the  displacement  gradient  with  respect  to  the 
coordinates  of  the  unstrained  system,  duk/ dr(°  or  V0u. 

If  the  one-electron  wavefunction  is  expanded  according  to  the  envelope  function  approximation 
in  unstrained  coordinates,  then 


V'aCxC r°))  =  */>a( r°)  =  . 


(235) 


Applying  £jkp*}pk  to  a  term  in  the  this  envelope  function  expansion  yields 

ejkPjPkffn'\r°^Um<r0)  =  £jk\PjPkUm'(.r° )  +  U m'  (X°)Pj  Pk\fffi  (r°) 

+  £jk{[PjUm'(r°)]Pk  +  [PkUm'(.r°)]Pj}ffi(.r°)  , 

=  ejklPjPkUm’( r°)  +  h2 Umf  (r0)kj  kk\fffi (r°) 

+  2hejk [pj  Um> (j°)]kkfffi (r°)  , 

where  kf  =  pf  /h.  Following  Burt  (57),  the  following  expansions  may  be  made: 

P]  U  m'  (r°)  =  YtmV^'lUm  (r°)  , 


and 


P'jPkUm' 


Um( r°)  ■ 


The  effective  potential  may  be  expanded  as  a  Taylor  series  about  €jk  —  0  and  o)jk  —  0. 

^ext,eff(x(r°).  Wsys) 

*  Wf(r°>{R}sys)  +  ^(r°'{R)sys)^  +  ^(r°,  {R}sys)<»i7  , 


where 


. .  _  ^  ^ext.eff 

iJ-^~ 


(jj  _  ^^ext.eff 


£ij= 0 

<wy= 0 


and  V-f  = 

1  daiij 


£ij  =  0 

0>ij= 0 


(236) 

(237) 

(238) 


(239) 


(240) 


If  the  system  is  subject  to  only  small  rotations,  then  a  stationary  test  charge  is  unlikely  to 
experience  much  of  a  change  in  potential  as  the  system  rotates  about  it,  and  Vff  is  likely  to  be 

negligible.  If  one  presumes  that  equation  172  still  holds  for  the  undefonned  configuration  (i.e., 
with  r°  in  place  of  r),  and  that 


Vtj{ r°,  {R}sys)  *  ^  0i(r° Wr°'  ’ 


(241) 


then  one  may  treat  Fj;(r°,  (R}sys)  much  as  Burt  (57)  treats  14xt,eff  (r°,{R}Sys)  in  equation  173; 
that  is, 
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Kj7(r°,{R}sys)^(r°)  « 

Zi  0,(r°)  Z„.  Em  [J„c  ul,(r°)V„(r°,  {R}L)ym(r»)d3r]  f%\ r»)  • 
Applying  the  envelope  function  approximation,  then,  leads  to  the  following  equations: 
2m^toth0.WSyS.«)/ir)(r0)  *  Elffl (r«)  , 

where 

’,(r°.{R}sys.e)  =  Z,  0,(r)/?E,  (»&,€)  , 


and 


ft 


/?^({R}L,€)  =  H^;«R}t,)  -  2e,*)pKft°  +  *liD% 


»0'fe) 

'm 


+  2m„  ^  2€J*)kjkk  2 


—  (ekj  +  a)kj) 


[h  po»>y 

m  / 


m 


(242) 


(243) 

(244) 


(245) 


andH*r;({R}L  e)  are  almost  analogous,  with  the  former  obtainable  from  the 
latter  by  making  the  replacement  k°  — >  k°  and  adding  a  gradient  tenn  that  is  the  last  term  in  the 
above  equation.  This  gradient  tenn  will  be  zero  if  pffim  is  zero,  which  happens  when  Um( r°) 
and  Um'( r°)  are  linear  combinations  of  X(r°),  T(r°),  and  Z( r°),  e.g.,  the  class  A  basis 
functions  of  the  Luttinger-Kohn  formulation,  or  when  Um(r°)  and  Um>  (r°)  are  multiples  of 
S(r°).  lfUm'(r° )  =  S(r°),  Um(r°)  =  Z(r°),  and  A:  =  3,  then  -  ih is  Kane’s 
parameter,  PK ane,  from  equation  130.  Usually,  the  gradient  term  is  neglected,  but  the  envelope 
function  formulation  by  Zhang  (101)  has  incorporated  it. 

Previously,  it  was  mentioned  that  the  envelope  function  approximation  may  not  be  valid  in 
homogeneously  strained  regions,  since  local  periodicity  may  be  effectively  lost  in  such  regions. 
However,  the  envelope  function  expansion  in  equation  235  is  in  tenns  of  the  undeformed 
coordinates,  and  this  corrects  for  loss  of  local  periodicity.  This  envelope  function  expansion  may 
be  rewritten  as 


Va  (0  ='Ydffi{x  1  (r)  )Um’(x  1  (r) )  .  (2 

m' 

If  the  original  undeformed  system  was  composed  of  regions  that  are  locally  periodic,  then  the 
periodicity  of  Umi  in  the  undeformed  coordinate  r°  reflects  this.  The  function  Um>  (x1  (r)), 
though,  is  not  in  general  periodic  in  r,  and  the  argument  X-1(r)  allows  this  function  to  track 
departures  from  local  periodicity.  This  is  illustrated  in  figure  27. 
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Figure  27.  Schematic  of  a  function  t/„(x  '(r)) 

superimposed  over  a  distorted  lattice  where 
the  spacing  between  neighboring  atoms  is  no 
longer  the  same.  The  inverse  deformation 
map  x  1  allows  the  function  to  track  the 
distortions  in  the  lattice. 


5.  Valence  Band  Offset 


If  one  is  determining  the  band  structure  of  a  single  material,  the  choice  of  zero-energy  datum  is 
largely  immaterial.  In  the  example  band  structures  shown  in  figures  4,  9,  and  18,  the  valence 
band  maximum  was  arbitrarily  chosen  to  be  this  zero-energy  datum.  However,  if  two  materials 
are  brought  together,  their  one-electron  energies  need  to  be  determined  with  respect  to  a  common 
energy  datum,  and  this  is  often  expressed  in  terms  of  the  valence  band  offset.  That  is,  if  the 
valence  band  offset  between  materials  A  and  B  is  EVB0,  then  the  valence  band  maximum  in  bulk 
material  B  is  set  to  be  higher  than  the  valence  band  maximum  of  A  by  £Vbo-  This  band  offset 
may  be  determined  by  a  variety  of  theoretical  and  experimental  means  (102),  and  once  it  is 
determined,  it  may  be  used  to  modify  the  parameters  for  the  various  methods.  For  example,  if 
the  fitting  parameters  for  each  material  lead  to  band  structures  where  the  valence  band  maximum 
for  each  material  is  zero,  then  in  the  tight  binding  method,  the  diagonal  elements  Homom  of  the 
Hamiltonian  matrix  for  the  atoms  in  material  B  are  replaced  by  Homom  +  EVB0.  Similarly,  for 
the  envelope  function  method,  diagonal  elements  R}m)  would  be  replaced  by 

H-mm  ({R}S>)  +  Fvbo.  This  only  has  to  be  done,  however,  if  the  fitting  parameters  do  not  already 
take  the  valence  band  offset  parameters  into  account.  For  example,  the  fitting  parameters  from 
Boykin  et  al.  (33)  for  GaAs  and  InAs  lead  to  valence  band  maxima  of  0.0  and  0.22  for  each 
material,  respectively,  so  an  additional  offset  should  not  be  added  to  Homom  when  these 
parameters  are  used. 


6.  Example  Results  and  Discussion 


With  the  help  of  available  data  in  the  literature,  we  next  present  results  of  comparisons  among 
the  empirical  methods  described  previously  for  model  systems  of  (1)  bulk  GaAs,  (2)  a  slab  of 
InAs,  (3)  an  AlAs/GaAs/AlAs  quantum  well,  and  (4)  an  InAs  quantum  dot  embedded  in  an 
InGaAs  matrix.  We  use  the  plane -wave  version  of  the  empirical  pseudopotential  method  of 
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section  3.1  and  the  Slater-Koster  tight-binding  method  of  section  3.2.  The  bulk  crystal  results 
pertain  to  the  effects  of  strain  on  the  electronic  structure  at  the  T  point.  These  results  are 
compared  with  the  corresponding  results  predicted  from  a  form  of  the  k  ■  p  method  that  accounts 
for  the  effects  of  strain  on  the  valence  bands,  as  discussed  in  sections  4.2  and  4.5,  combined  with 
an  effective  mass  method  for  bulk  crystals,  discussed  in  section  4.4,  to  account  for  the 
conduction  band.  The  particular  strain  Hamiltonian  used  for  the  valence  bands  is  equation  222, 
while  the  change  in  the  conduction  band  minimum  with  strain  is  accounted  through  equation 
227.  The  InAs  slab  is  modeled  through  the  tight-binding  method  under  arbitrary  homogeneous 
strains.  And  the  strain-free  GaAs  quantum,  well  clad  by  two  AlAs  layers,  is  modeled  by  both  the 
tight-binding  and  envelope-function  methods.  Finally,  results  from  the  tight-binding  and 
envelope-function  methods  are  shown  for  an  InAs  quantum  dot  embedded  in  a  layer  of 
Ino.4Gao.6As. 

Figure  28  shows  a  comparison  of  results  from  the  empirical  pseudopotential  method  (EPM)  by 
Mader  and  Zunger  ( 47)  and  a  k  ■  p  fonnulation  by  Van  de  Walle  (103),  both  applied  to  a  bulk 
crystal  of  GaAs.  The  k  ■  p  formulation  used  here  is  equivalent  to  using  the  Chuang  k  ■  p 
Hamiltonian  in  equation  161  with  the  deformation  potential  values  from  Van  de  Walle,  together 
with  an  effective  mass  approximation  for  the  conduction  band.  The  first  subfigure  shows  the 
energies  of  conduction  and  valence  electrons  for  wave  vector  k  =  0,  that  is,  the  band  edges 
calculated  by  the  two  methods,  while  the  second  subfigure  shows  the  band  gaps  determined  from 
those  energies.  Agreement  between  the  two  methods  is  far  better  for  the  band  gap  results  than 
for  the  electronic  energies  themselves,  which  illustrates  the  problem  discussed  by  Williamson 
and  Zunger  (41),  who  found  that  without  an  explicitly  strain-dependent  pseudopotential,  the 
changes  in  band  gaps  could  be  fit  to  experiment  but  not  the  changes  in  the  band  edges 
themselves.  The  EPM  formulation  of  Mader  and  Zunger  does  not  contain  any  explicitly  strain- 
dependent  terms,  unlike  the  formulation  of  Kim  et  al.  (50).  The  variation  with  strain  of 
electronic  energies  due  to  this  latter  formulation  is  shown  in  figure  29,  along  with  the 
corresponding  k  ■  p  results  of  Van  de  Walle.  The  latter  formulation  better  captures  the  changes 
in  electronic  energies  due  to  hydrostatic  strain. 
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Figure  28.  For  bulk  GaAs,  the  variation  with  hydrostatic  strain  of  (a)  the  energies  of  conduction  and  valence 
electrons  for  wave  vector  k  =  0  and  (b)  the  band  gap,  i.e.,  the  difference  between  the  energies  of 
conduction  and  valence  electrons,  for  wave  vector  k  =  0.  The  electronic  energies  are  calculated 
using  the  empirical  pseudopotential  method  (EPM)  as  formulated  by  Mader  and  Zunger  (47)  and 
using  a  k  ■  p  formulation  from  Van  de  Walle  (103). 


Figure  29.  For  bulk  GaAs,  the  variation  with  hydrostatic 
strain  of  the  energies  of  conduction  and 
valence  electrons  for  wave  vector  k  =  0, 
where  the  electronic  energies  are  calculated 
using  the  EPM  as  formulated  by  Kim  et  al. 
(50)  and  using  a  k  ■  p  formulation  from  Van 
de  Walle  (705). 
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As  shown  in  figure  30,  the  strain-dependent  EPM  formulation  of  Kim,  Wang,  and  Zunger  also 
tracks  changes  in  electronic  energies  due  to  strain  states  besides  hydrostatic  strain,  such  as  strain 
due  to  a  biaxial  stress  state  in  the  plane  nonnal  to  the  (001)  crystal  direction  of  GaAs.  Here,  the 
z-direction  is  taken  along  (001),  and  the  nonnal  stresses  point  along  the  x-  and  y-axes  and  are 
equal  in  magnitude.  These  stresses  are  either  both  tensile  or  both  compressive,  so  that 
exx  —  eyy-  and  the  value  of  ezz  is  detennined  from  the  Poisson  effect.  The  first  subfigure 
compares  the  EPM  results  with  a  spin-free  k  ■  p  formulation  that  accounts  for  the  heavy  and 
light-hole  bands  (81),  along  with  the  effective  mass  approximation  for  the  conduction  band  used 
earlier.  The  second  subfigure  uses  Van  de  Walle’s  k  ■  p  formulation,  which  includes  spin 
effects.  This  particular  strain-dependent  EPM  fonnulation  does  not  account  for  spin-orbit 
coupling,  and  so  it  agrees  better  with  the  spin-free  k  ■  p  formulation,  especially  with  regard  to 
the  light-hole  band.  Without  spin-orbit  coupling,  the  variation  of  the  light-hole  band  edge  with 
strain  is  linear.  The  tight-binding  fonnulation  of  Boykin  et  al.  (33) — which  does  account  for 
spin — does  capture  the  nonlinearity  in  the  variation  from  strain  due  to  spin-orbit  coupling. 


(a)  (b) 


Figure  30.  For  bulk  GaAs,  the  variation  with  biaxial  stress-induced  strain  of  the  energies  of  conduction  and 

valence  electrons  for  wave  vector  k  =  0,  where  the  electronic  energies  are  calculated  with  the  EPM  as 
formulated  by  Kim  et  al.  (50)  and  also  (a)  a  spin-free  k  ■  p  formulation  and  (b)  a  k  ■  p  formulation 
with  spin.  Both  k  ■  p  formulations  use  the  deformation  potential  values  from  Van  de  Walle  (103). 

The  strain  along  the  horizontal  axis  is  the  strain  in  the  plane  normal  to  the  (001)  crystal  direction. 
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Figure  31.  For  bulk  GaAs,  the  variation  with  biaxial 
stress-induced  strain  of  the  energies  of 
conduction  and  valence  electrons  for  wave 
vector  k  =  0,  where  the  electronic  energies 
are  calculated  with  the  tight-binding 
formulation  of  Boykin  et  al.  (Ji)  and  also  a 
k  ■  p  formulation  from  Van  de  Walle  (103). 

The  strain  along  the  horizontal  axis  is  the 
strain  in  the  plane  normal  to  the  (00 1 )  crystal 
direction. 

To  further  illustrate  the  effect  of  strain  on  electronic  structure  results,  the  band  structures  of  a 
slab  of  InAs  about  1.2  nm  thick  subjected  to  in-plane  stresses  are  shown  in  figure  32.  The 
stresses  are  along  the  (100)  and  (010)  directions,  leading  to  a  strain  of  f||  along  these  same 
directions.  No  stress  is  imposed  along  the  (001)  direction.  One-particle  energies  are  shown  for 
electrons  with  Bloch  wave  vectors  pointing  along  the  (100)  and  (110)  crystal  directions.  The 
band  structure  is  determined  through  the  tight-binding  method.  Dangling  bonds  are  again 
tenninated  with  hydrogens.  Essentially,  this  is  an  idealized  quantum  well  where  electrons  are 
confined  not  by  semiconductor  layers  but  by  vacuum,  and  the  strain  is  allowed  to  vary  arbitrarily 
rather  than  be  fixed  by  the  lattice  mismatch  between  the  material  of  the  well  and  its  surrounding 
semiconductor.  The  band  gap  narrows  with  increasing  in-plane  tensile  strain  and  widens  with 
increasing  in-plane  compressive  strain,  and  increasing  in-plane  tensile  strain  also  leads  to  the 
valence  bands  bunching  together  and  crossing  over  one  another. 
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Figure  32.  Band  structure  of  a  1.2-nm  thick  InAs  slab  subject  to  stresses  along  the  (100)  and  (010)  directions, 
leading  to  a  strain  of  eN  along  these  same  directions.  No  stress  is  imposed  along  the  (001)  direction. 

The  &||  values  to  the  left  of  zero  are  the  negative  magnitudes  of  Bloch  wavevectors  pointing  along 
crystal  direction  (100),  while  the  lq  values  to  the  right  of  zero  are  the  magnitudes  of  Bloch 
wavevectors  pointing  along  crystal  direction  (110). 

A  comparison  of  the  results  of  the  tight-binding  and  envelope  function  methods  is  shown  in 
figure  33.  These  methods  were  used  to  find  the  one-particle  energies  of  electrons  with  Bloch 
wave  vectors  pointing  along  the  (100)  and  (110)  crystal  directions  in  a  5 -nm- thick  GaAs 
quantum  well  sandwiched  between  30-nm-thick  layers  of  AlAs.  Two  sets  of  tight-binding 
results  were  generated  with  NEM05  modeling  code  (35),  one  using  parameters  from  Jancu  et  al. 
(104)  and  one  using  parameters  from  the  NEM05  material  database.  The  GaAs  parameters  from 
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Valence  band  dispersion,  AlAs/GaAs/AIAs  QW 


Valence  bands,  EFA  -  TB,  Jancu  params  +  TB,  Klimeck  params 

(a) 


Conduction  band  dispersion,  AlAs/GaAs/AIAs  QW 


Conduction  bands,  EFA  -  TB,  Klimeck  params  x 

TB,  Jancu  params  + 

(b) 


Figure  33.  Comparison  of  results  from  tight-binding  and  envelope  function  methods  for  a  5-nm-thick  GaAs 

quantum  well  sandwiched  between  30-nm-thick  layers  of  AlAs.  Valence  bands  are  shown  in  (a)  and 
conduction  bands  in  (b).  Two  sets  of  tight-binding  parameters  are  used,  one  from  Jancu  et  al.  (104)  and 
one  from  used  in  the  material  database  of  the  tight-binding  software  code  NEM05  and  attributed  to  one 
of  its  authors,  G.  Klimeck  (35).  The  8*8  k  ■  p  Hamiltonian  used  by  the  envelope  function  calculations 
may  be  found  in  the  documentation  of  nextnano  (105)  or  Andlauer  (106).  The  k\\  values  to  the  left  of 
zero  are  the  negative  magnitudes  of  Bloch  wavevectors  pointing  along  crystal  direction  (100),  while  the 
&n  values  to  the  right  of  zero  are  the  magnitudes  of  Bloch  wavevectors  pointing  along  crystal  direction 
(110). 

this  database  may  be  found  in  Boykin  et  al.  (33),  and  the  AlAs  parameters  from  this  database  are 
shown  in  table  4.  A  valence  band  offset  of  between  AlAs  and  GaAs  was  needed  when  the 
parameters  from  Jancu  et  al.  were  used,  and  the  valence  band  maximum  of  AlAs  was  taken  to  be 
0.5  eV  lower  than  that  of  GaAs.  However,  the  parameters  from  the  NEM05  database  already 
incorporated  a  valence  band  offset.  Boundary  conditions  were  periodic  only  along  the  directions 
in  the  plane  of  the  well,  (100)  and  (010),  and  dangling  bonds  were  tenninated  with  hydrogen 
atoms.  The  envelope  function  results  were  generated  from  nextnano  simulation  software  using 
an  8  x  8  k  ■  p  Hamiltonian  matrix  documented  on  the  nextnano  web  site  (105)  and  by  Andlauer 
(106).  Agreement  between  the  two  methods  is  good  for  the  valence  bands,  especially  the  top  two 
bands.  For  the  conduction  bands,  the  tight-binding  results  using  the  parameters  of  Jancu  et  al. 
only  agree  with  the  envelope  function  results  for  small  wave  vector  values,  and  the  tight-binding 
results  using  the  other  set  of  parameters  lead  to  a  conduction  band  minimum  about  0. 1  eV  less 
than  that  predicted  from  the  envelope  function  method.  For  the  valence  bands,  a  similar 
agreement  between  tight-binding  and  envelope  function  results  has  been  demonstrated  by  de 
Franceschi  et  al.  (107),  who  compared  tight-binding  results  using  the  parameters  of  Jancu  et  al. 
with  results  from  envelope  function  results  using  a  6  x  6  k  ■  p  Hamiltonian  matrix. 
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Table  4.  Tight-binding  parameters  for  AlAs  from  the  material 

database  of  NEM05  (35).  These  parameters  are  values  of 
E®}  and  V (,°') ,  from  equations  93  and  97,  and  the 

o^m  o  m  ,om,A  n  ’ 

spin  parameter  Am  in  equation  88  for  o' ,  o  £  (s,  p,  d,  5*} 
and  m',  m  £  (Al,  As}. 


p(0) 

csZA1 

0.79695 

1/  (°)  y(°) 

sAl,sAs,er  ’  sAs,sAl,er 

-1.64584 

p  (0) 

CsZas 

-5.17012 

l/(°)  1/1°) 

pAl,pAs,er’  pAs,pAl,er 

4.53156 

p(°) 

PZA  1 

6.63291 

T/(°)  y(°) 

pAl,pAS,7T’  pAS,pAl,7T 

-1.86816 

p(°) 

cpzAs 

4.39708 

^(°)  y® 

^dAl.dAs.er’  ^dAs,dAl,er 

-1.97058 

p 1°) 
CdZA1 

12.92120 

y(°)  y® 

^dAl,dAs,7r’  ^dAs,dAl,7r 

1.67733 

p(°) 

CdZAs 

13.13880 

y(°)  y® 

^dAl,dAs,<5’  vdAs,dA\,S 

-1.58868 

p(°) 

VzA, 

24.16590 

y(°)  y(°) 

^s*Al,s*As,cr’  ^s*As,s*Al,cr 

-2.84245 

p(°) 

C‘s*ZA, 

19.80470 

— 

— 

y(°) 

sAl,pAs,cr 

2.95309 

y(°) 

sAs,pAl,cr 

3.02223 

y(°) 

GAl,dAs,er 

-2.64111 

y(°) 

^sAs,dAl,er 

-3.03196 

y(°) 

^sAl,s*As,cr 

-1.88341 

y(°) 

GAs,s*Al,ff 

-2.78690 

r/(°) 

^pAl,dAs,er 

-1.02836 

y(°) 

^pAs,dAl,er 

-2.47345 

y(°) 

pAl,dAs,7r 

2.86419 

y(°) 

pAs,dAl,7r 

2.52741 

y(°) 

vs*Al,pAs,a 

1.30469 

y(°) 

^s*As,pAl,er 

1.92174 

y(°) 

vs*  Al, dAs,ff 

-1.73510 

y(°) 

^s*As,dAl,er 

-1.84300 

^A1 

0.01586 

'Us 

0.17386 

A  comparison  of  the  results  of  the  tight-binding  and  envelope  function  methods  is  shown  in 
figure  33.  These  methods  were  used  to  find  the  one-particle  energies  of  electrons  with  Bloch 
wave  vectors  pointing  along  the  (100)  and  (110)  crystal  directions  in  a  5 -nm- thick  GaAs 
quantum  well  sandwiched  between  30  nm-thick  layers  of  AlAs.  Two  sets  of  tight-binding  results 
were  generated  with  NEM05  (35),  one  using  parameters  from  Jancu  et  al.  (104)  and  one  using 
parameters  from  the  NEM05  material  database.  The  GaAs  parameters  from  this  database  may 
be  found  in  Boykin  et  al.  (33),  and  the  AlAs  parameters  from  this  database  are  shown  in  table  4. 
A  valence  band  offset  of  between  AlAs  and  GaAs  was  needed  when  the  parameters  from  Jancu 
et  al.  were  used,  and  the  valence  band  maximum  of  AlAs  was  taken  to  be  0.5  eV  lower  than  that 
of  GaAs.  However,  the  parameters  from  the  NEM05  database  already  incorporated  a  valence 
band  offset.  Boundary  conditions  were  periodic  only  along  the  directions  in  the  plane  of  the 
well,  (100)  and  (010),  and  dangling  bonds  were  terminated  with  hydrogen  atoms.  The  envelope 
function  results  were  generated  from  nextnano  using  an  8  x  8  k  ■  p  Hamiltonian  matrix 
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documented  on  the  nextnano  web  site  (105)  and  by  Andlauer  (106).  Agreement  between  the  two 
methods  is  good  for  the  valence  bands,  especially  the  top  two  bands.  For  the  conduction  bands, 
the  tight-binding  results  using  the  parameters  of  Jancu  et  al.  only  agree  with  the  envelope 
function  results  for  small  wave  vector  values,  and  the  tight-binding  results  using  the  other  set  of 
parameters  lead  to  a  conduction  band  minimum  about  0. 1  eV  less  than  that  predicted  from  the 
envelope  function  method.  For  the  valence  bands,  a  similar  agreement  between  tight-binding 
and  envelope  function  results  has  been  demonstrated  by  de  Franceschi  et  al.  (107),  who 
compared  tight-binding  results  using  the  parameters  of  Jancu  et  al.  with  results  from  envelope 
function  results  using  a  6  x  6  k  ■  p  Hamiltonian  matrix. 

Comparisons  between  tight-binding  and  envelope  function  results  have  also  been  made  for  a 
dome-shaped  InAs  quantum  dot  by  Sengupta  et  al.  (108).  The  dot  was  5  nm  in  height,  20  nm  in 
diameter,  and  embedded  in  a  5-nm  layer  of  Ino.4Gao.6As,  which  in  turn  was  sandwiched  between 
30-nm  layers  of  GaAs.  The  strain  in  the  quantum  dot  was  determined  through  two  different 
atomistic  valence  force  field  (VFF)  models,  the  harmonic  Keating  model  (1 09)  and  an 
anharmonic  model  (110).  Results  from  the  tight-binding  and  envelope  function  methods  are 
shown  in  table  5.  The  theoretical  results  were  compared  to  experimental  results  for  a  similar 
quantum  dot  whose  peak  photoluminescence  occurred  at  the  wavelength  A  —  1.52  jum  (111), 
corresponding  to  a  band  gap  of  about  Eg  «  2nhc/A  ~  0.82  eV.  Agreement  between  the  two 
theoretical  methods,  tight-binding  and  the  envelope  function  approximation,  is  good  for  both 
methods  of  determining  the  strain,  but  agreement  with  experiment  is  much  better  when  the 
anharmonic  VFF  model  is  used  to  determine  strain. 

Table  5.  Results  for  the  quantum  dot  studied  by  Sengupta  et  al. 

(108).  Dot  is  5  nm  in  height,  20  nm  in  diameter,  and 
embedded  in  a  5-nm  layer  of  In0.4Ga0.6As,  sandwiched 
between  two  30-nm  GaAs  layers. 


Model 

Band  Gap  (eV) 

Envelope  function,  harmonic  VFF 

1.063 

Tight-binding,  harmonic  VFF 

1.040 

Envelope  function,  anharmonic  VFF 

0.885 

Tight-binding,  anharmonic  VFF 

0.828 

7.  Conclusion 


Various  methods  for  estimating  electronic  structure  have  been  discussed,  along  with  discussions 
of  how  strain  is  incorporated  into  these  methods.  We  did  not  find — among  the  methods  both 
included  and  not  included  in  this  survey — a  generalized  and  explicit  treatment  of  continuum 
deformability.  Thus  a  gap  exists  between  the  methods  suited  to  handle  slowly  varying  elastic 
fields  and  fully  atomistic  approaches.  Atomistic  approaches  offer  the  greatest  generality  but  lack 
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an  ability  to  treat  deformation  explicitly  without  recourse  to  a  minimization  of  atomic  scale 
forces  and  its  attendant  computational  expense.  They  categorically  require  a  full  recalculation 
over  all  atomistic  degrees  of  freedom  based  on  each  new  configuration  of  nuclei,  which  therefore 
scales  according  to  the  number  of  atoms,  electrons,  or  electronic  wavefunctions  regardless  of  the 
degree  to  which  the  deforming  material  exhibits  continuum-like  behavior.  Continuum-based 
methods  incorporate  notions  of  stress,  but  rely  exclusively  on  linear  elasticity  or  infinitesimal 
strain  theories.  Elastic  fields  are,  at  best,  mildly  inhomogeneous.  Namely,  finite  deformation 
effects  such  as  near  dislocation  or  vacancy  cores,  particularly  for  complex  lattices,  are  not 
presently  possible.  Results  for  finite  strain  effects,  therefore,  require  a  separate  atomistic 
calculation  for  verification — first  to  minimize  the  energy  of  the  nuclear  configuration  in  the 
deformed  state  and  secondly,  in  some  cases,  to  detennine  the  resulting  electronic  structure.  For 
systems  whose  sizes  simultaneously  need  to  be  sufficiently  large  to  capture  the  convergence  of 
the  elastic  field  solution  to  the  bulk  limit  while  also  providing  the  required  resolution  near 
features  of  interest,  such  as  at  interfaces  and  defects,  single  point  calculations  alone  may  be 
computationally  costly  and  calculations  over  multiple  mechanically  deformed  states  may  be 
prohibitive. 

In  ab  initio  methods,  strain  is  taken  into  account  simply  through  the  positions  of  the  atoms  in  the 
system  being  simulated.  However,  in  atomistic  empirical  methods,  the  positions  of  the  atoms  are 
not  always  sufficient  to  fully  take  the  effects  of  strain  into  account.  In  the  empirical 
pseudopotential  method,  a  strain-dependent  prefactor  in  the  atomic  form  factors  has  been  used  to 
better  account  for  the  effects  of  strain  on  band  edges  (41).  In  the  tight-binding  method,  atomic 
positions  are  usually  sufficient  to  account  for  the  effects  of  strain.  However,  if  there  is  long- 
range  charge  transfer,  such  as  that  due  to  piezoelectric  material  properties,  then  the  resulting 
electric  field  has  to  be  treated  like  an  external  field  in  these  atomistic  methods,  and  in  the  case  of 
piezoelectricity,  this  field  is  determined  through  continuum  mechanical  calculations  (34,  78). 

The  k  ■  p  and  envelope  function  methods  are  empirical  electronic  structure  methods  that  do  not 
depend  at  all  on  the  positions  of  the  atoms  in  the  simulation.  Rather,  the  potential  felt  by  an 
electron  in  a  region  of  the  system  composed  of  some  material  is  assumed  to  be  the  potential  in  a 
bulk  crystal  of  the  material.  Because  of  this,  they  are  less  computationally  expensive  than  the 
atomistic  methods,  but  the  positions  of  the  atoms  can  no  longer  be  used  to  take  strain  into 
account.  Instead,  a  strain-dependent  term  is  added  to  the  one-electron  Hamiltonian.  This  term 
was  originally  derived  by  Bir  and  Pikus  (I)  for  an  infinite  bulk  subjected  to  homogeneous  strain. 
However,  the  tenn  can  be  derived  in  an  alternative  fashion  for  more  general  states  of  strain, 
though  the  gradient  of  the  strain  still  needs  to  be  small  (100). 

Finally,  examples  of  the  results  from  the  empirical  electronic  structure  methods  discussed  have 
been  shown.  Generally,  when  the  methods  are  applied  to  the  same  problem  and  each  is 
examined  within  the  assumptions  unique  to  each  method,  the  results  are  mostly  in  agreement. 
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Appendix  A.  Finding  Reciprocal  Lattice  Vectors  Within  a  Given  Energy 
Cutoff 


where  Gcut  is  the  magnitude  of  the  largest  reciprocal  lattice  vector  satisfying  the  cutoff  criterion. 
Since  G  =  Y,f=i  ni  b*,  the  largest  possible  value  of  rij  is1 


j-max  Gcut 

i  |b;|  ’ 


(A-2) 


A  method  for  finding  the  reciprocal  lattice  vectors  satisfying  the  energy  cutoff  criterion  for 
k  =  0  is  shown  in  figure  A-l.  Once  these  vectors  have  been  found,  another  method,  shown  in 
figure  A-2,  may  be  used  to  find  reciprocal  lattice  vectors  satisfying  the  energy  cutoff  criterion  for 
nonzero  k.2 


'  Kohanoff,  J.  Electronic  Structure  Calculations  for  Solids  and  Molecules',  Cambridge  University  Press:  New  York,  2006. 
2 

“Varga,  K.;  Driscoll,  J.  A.  Computational  Nanoscience',  Cambridge  University  Press:  New  York,  2011. 
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£cu.  is  given. 

Primitive  reciprocal  lattice  vectors  b, 
are  given. 


Gcut  =  \/2meEcM/h2 

for  i  =  1  to  3  do 
«,max  =  Gcut/|b,| 

end  for 

Initialize  list  GvecList. 

for  «|  =  -«'jnax  to  n'jnax  do 
for  ni  =  —  n™ax  to  «™ax  do 
for  =  — «”lax  to  n”,ax  do 

G  =  « i  b  i  +ri2bi  +«3b3 
if  |G|  <  Gcut  then 
Add  G  to  GvecList. 

end  if 

end  for 
end  for 

end  for 


Figure  A-l.  Pseudocode  for  finding 

reciprocal  lattice  vectors  G  that 
satisfy  the  energy  cutoff  criterion 
when  k  =  0. 
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EcM  is  given. 

Wavevectors  k| .  ki, . . . ,  k,yk  are  given. 


Also  given  is  the  list  GvecList  of  reciprocal  lattice  vectors  whose 
magnitudes  are  less  than  Gcut  =  \]2meEcui/tr . 


for  k  in  k|  .ki, _ kA/k  do 


Initialize  list  NewGvecList{ k). 
for  G  in  GvecList  do 


if  |k  +  G|  <  \j2meEcM/tr  then 
Add  G  to  N ewGvecList  ( k ) . 

end  if 


end  for 
end  for 


Figure  A-2.  Pseudocode  for  finding  reciprocal  lattice  vectors  G  that 
satisfy  the  energy  cutoff  criterion  for  nonzero  k. 
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Intentionally  left  blank. 
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Appendix  B.  Calculating  Strain  From  Atomic  Positions 


Figure  B-l.  Before  the  strain  is  applied,  the 

atoms  are  taken  to  be  in  their  ideal 
tetrahedral  positions,  with  the 
lighter  colored  atom  at  the  center 
of  the  tetrahedron.  After  the  strain, 
the  tetrahedron  is  distorted,  and 
vectors  T°,  U°,  and  V°  become  T, 

U,  and  V. 


There  is  no  unique  method  to  detennine  the  strain  from  atomic  coordinates.  Saito  and  Arakawa1 
and  Steiger  et  al.2  use  the  method  of  Pryor  et  al.3 * *  to  determine  the  strain.  Figure  B-l  illustrates  the 
positions  of  some  atoms  before  and  after  strain.  The  atoms  are  tetrahedrally  coordinated,  with  the 
atom  shown  in  lighter  color  being  at  the  center  of  the  tetrahedron.  After  the  strain  is  applied,  this 
tetrahedron  is  distorted.  The  vectors  T°,  U°,  and  V°  are  aligned  to  the  edges  of  the  ideal  tetrahedron, 
while  the  vectors  T,  U,  and  V  are  aligned  to  the  edges  of  the  distorted  one.  The  relationships 
between  these  sets  of  vectors  and  the  strain  are  taken  to  satisfy  the  matrix  equation 
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which  can  be  inverted  to  find  the  strain,  so  that 
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(B-l) 


(B-2) 


Equation  B-l  is  equivalent  to  the  following  set  of  tensor  equations: 


'Saito,  T.;  Arakawa,  Y.  Electronic  Structure  of  Piezoelectric  In(0.2)Ga(0.8)N  Quantum  Dots  in  GaN  Calculated  Using  a 
Tight-Binding  Method.  PhysicaE:  Low-Dimensional  Systems  and  Nanostructures  2002,  15,  169-181. 

J 

Steiger,  S.;  Povolotskyi,  M.;  Park,  H.-H.;  Kubis  T.;  Klimeck,  G.  NEM05:  Parallel  Multiscale  Nanoelectronics  Modeling 
Tool.  IEEE  Transactions  on  Nanotechnology  2011,  10,  1464-1474. 

3 

Pryor,  C.;  Kim,  J.;  Wang,  L.  W.;  Williamson,  A.  J.;  Zunger,  A.  Comparison  of  Two  Methods  for  Describing  the  Strain 

Profiles  In  Quantum  Dots.  J.  of  Applied  Physics  1998,  S3,  2548-2554. 


99 


Ti  =  (Sij  +  6lj)t;\ 

ut  =  (Sij  +  e^Uj0, 


and 


v,  =  (S,j  +  e,j)vf  •  (B-3) 

Here,  a  repeated  index  implies  summation  from  1  to  3  over  that  index.  These  equations  are  a 
more  specific  version  of  the  Cauchy-Born  rule,4  where  5^  +  has  replaced  the  deformation 
gradient  Fi;-.  They  can  be  interpreted  as  a  version  of  the  Cauchy-Born  rule  that  neglects  rotation, 


since 


Fij  -  sij  +  §r<7  -  S^  +  (eij  +  <x>ij)  »  Sij  +  £ij , 


dr 


(B-4) 


where  iq  is  the  displacement,  r®  is  the  position  vector  expressed  in  the  coordinates  of  the 
undeformed  system,  and  et j  and  o>LJ  are 


eij=~[ 


lfdiii  duj N 
2  \dr°  +  dr?  y 


and 


1/dUi  duj'' 

IJ  2  1  dr?  dr? , 


(B-5) 


j  i 

Since  oi>i  .■  is  neglected,  it  is  implicitly  assumed  that  Fu  is  symmetric.  If  it  is  not,  the  values  of  e 


ij 


determined  from  method  of  Pryor  et  air  will  not  be  either. 


Another  method  to  determine  strain  from  changes  in  atomic  coordinates  treats  deviations  from 
the  Cauchy-Born  rule  as  residuals  to  be  minimized.  Horstemeyer  and  Baskes6  take  the 
minimizing  function  at  atom  i  to  be 

N 

0i(F)  =  ^]|Ry-FR?/|^O'),  (B-6) 

7  =  1 

where  W(j)  is  a  weighting  function,  F  is  the  deformation  gradient,  and  R°-  and  Rty  are  the 
vectors  connecting  atom  i  to  each  of  its  N  neighbors  in  the  undeformed  and  deformed 
configurations,  as  shown  in  figure  B-2.  If  the  Cauchy-Bom  rule  holds  exactly  at  atom  i,  then 
R ij  —  FR°;-  and  0;(F)  =  0.  \fW(J)  —  1,  the  value  of  F  that  minimizes  (pt  is7’8 


4Steinmann,  P.;  Elizondo,  A.;  Sunyk,  R.  Studies  of  Validity  of  the  Cauchy-Born  Rule  By  Direct  Comparison  of  Continuum 
and  Atomistic  Modelling.  Modelling  and  Simuation  in  Materials  Science  and  Engineering  2007,  15,  S271-S281. 

5Mase,  G.  T.;  Mase,  G.  E.  Continuum  Mechanics  for  Engineers’,  2nd  ed.,  CRC  Press:  Boca  Raton,  FL,  1999. 

Horstemeyer,  M.  F.;  Baskes,  M.  I.  Strain  Tensors  at  the  Atomic  Scale;  In  Multiscale  Phenomena  in  Materials — Experiments 
and  Modeling’,  Materials  Research  Society:  Boston,  2000. 

7Zimmermann,  J.  A.  Continuum  and  Atomistic  Modeling  of  Dislocation  Nucleation  at  Crystal  Surface  Ledges.  Ph.D.  Thesis, 
Stanford  University,  Stanford,  CA,  2000. 

g 

Zimmermann,  J.  A.;  Bammann,  D.  J.;  Gao,  H.  Deformation  Gradients  for  Continuum  Mechanical  Analysis  of  Atomistic 
Simulations.  International  J.  of  Solids  and  Structures  2009,  46,  238-253. 
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where 


(B-7) 


F  =  WjW,-1  , 


and 


N 

Wi  =^(Ri;®R?y), 

7  =  1 

N 

w2  =  £(R?;  ®  R?;)  ■ 

7  =  1 


Operator  0  is  the  dyad  product.  From  F,  several  different  strain  tensors  can  be  obtained. 


Figure  B-2.  Relative  positions  of  the  neighbors  of  atom  i  before  and 
after  the  strain  is  applied. 


(B-8) 
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